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Methods  of  investigating  internal  waves  from  bathythermograph 
information  are  discussed.  A  series  of  651  hourly  observations 
taken  at  Ocean  Station  "P"  is  analyzed  by  two  methods  to  determine 
the  characteristics  of  the  internal  waves  present.  One  method 
involves  a  study  of  the  fluctuations  in  depth  of  six  selected 
isotherms;  the  other  involves  a  study  of  the  fluctuations  in 
depth  of  the  top,  center,  and  bottom  of  the  thermocline.  An 
objective  method  is  used  to  determine  the  depths  of  these 
features  from  the  bathythermograph  data. 

Spectra  of  these  features  are  calculated  by  Fast  Fourier 
transform,  hanned,  and  smoothed,  and  these  spectra  are  compared 
with  previously  published  results  obtained  from  a  different 
method  of  analysis  of  the  same  time  series. 
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I.   INTRODUCTION 

Tidal  period  internal  waves  in  the  deep  ocean  were  observed 
as  long  ago  as  1907,  by  the  Swedish  oceanographer 
Pettersson  1909  .   Little  work  was  done  in  describing  these 
waves  or  investigating  their  cause  before  1950,  although  Stokes 
developed  the  theoretical  treatment  which  could  describe  internal 
waves --neglecting  geostrophic  effects--  in  1847.   This  theory 
was  used  by  Ekman  1904  in  explaining  the  "dead  water" 
phenomenon  in  Norwegian  fjords.   Although  internal  waves  have 
been  observed  to  exist  for  more  than  half  a  century,  only 
recently  have  mathematical  tools  and  equipment  been  developed 
which  enable  investigators  to  more  closely  analyze  these  waves  in 
terms  of  periodic  components  and  energy  distribution.   Although 
oceanographers  are  now  in  a  position  to  say  a  great  deal  about 
the  various  frequencies  which  are  present  in  a  spectrum  of 
internal  waves,  the  presence  of  many  of  these  frequencies 
and  the  causes  of  them  have  not  yet  been  explained. 

The  most  useful  tools  for  analysis  of  internal  wave  data 
are  Fourier  analysis  and  power  spectrum  analysis.   Fourier 
analysis,  although  a  tool  of  mathematicians  for  decades, 
required  great  expenditure  of  time  to  be  carried  out  before  the 
advent  of  high-speed  digital  computers,  and  has  only  recently 
been  brought  into  use  in  analyzing  internal  wave  records.   Power 
spectrum  analysis,  developed  and  refined  by  Blackman  and 
Tukey  |l959  has  been  used  by  electronic  and  communications 
engineers  for  several  years,  but  has  only  recently  been  applied 


to  ocean  wave  analysis.   Programming  either  of  these  methods  for 
use  on  a  digital  computer  makes  possible  analysis  of  large  amounts 
of  data  accurately  in  a  reasonably  short  time. 

In  order  to  describe  internal  waves,  it  is  at  least  necessary 
to  obtain  a  time  or  space  series  of  observations  for  analysis. 
The  most  common  method  of  obtaining  a  time  series  is  to  observe 
temporal  changes  in  thermal  structure  at  a  fixed  location.   This 
is  sometimes  done  from  weather  observation  ships  located  at  deep 
ocean  stations.  One  of  these,  Ocean  Station  "PTt,  located  in  the 
northeastern  Pacific  Ocean  at  145  West  Longitude,  50  North 
Latitude  and  maintained  by  the  Department  of  Transport  of  the 
Dominion  of  Canada,  has  been  making  bathythermograph  observations 
at  a  maximum  interval  of  12  hours  since  July  1952.  From  time  to 
time,  bathythermograph  observations  are  made  from  the  ship 
occupying  Station  "P"  at  four  hour  intervals,  or  more  frequently. 
The  resulting  bathythermograph  data  provides  an  excellent  time 
series  for  analysis  of  internal  waves. 

The  data  for  use  in  this  thesis  consisted  of  hourly 
bathythermograph  observations,  obtained  with  mechanical  bathy- 
thermographs from  2000  GMT,  5  August  1961  to  2200  GMT, 
1  September  1961,  a  total  of  651  hourly  observations.   The  data 
were  furnished  in  digital  format  on  punched  cards  by  the 
National  Oceanographic  Data  Center,  Washington,  D.  C.   Digital 
analysis  was  performed  on  an  IBM  360  computer  at  the  Naval 
Postgraduate  School  Computer  Facility,  Monterey,  California. 


B a thy thermograms  of  the  data  used  in  this  thesis  may  be  found 
in  Manuscript  Report  Series  No.  106,  published  by  the  Fisheries 
Research  Board  of  Canada,  Programmed  by  the  Canadian  Committee  on 
Oceanography. 


II.   OBTAINING  DATA  FOR  POWER  SPECTRUM  ANALYSIS 

One  of  the  most  commonly  used  methods  of  obtaining  data  for 
analysis  in  describing  internal  waves  has  been  following  the 
depth  fluctuations  of  selected  isotherms.   This  method  was  used 
by  Tabata  1965  at  Ocean  Station  "Pn  in  the  northeastern 
Pacific  Ocean,  and  by  Seiwell  1965   in  the  North  Atlantic  Ocean. 
Although  this  method  has  been  quite  successful  in  the  past,  its 
validity  depends  on  rather  specific  conditions.   Internal  waves 
propagate  along  isopycnal  surfaces,  and  have  their  greatest 
amplitude  at  a  density  discontinuity.   An  isothermal  surface 
will  closely  approximate  an  isopycnal  surface  for  propagation 
of  an  internal  wave  only  if  the  water  column  is  nearly  isohaline 
and  the  wave  is  of  small  amplitude .   In  this  case  one  can  approx- 
imate the  change  in  density  due  to  salinity  and  the  change  in 
density  due  to  pressure  fluctuations  by  zero.   In  choosing  the 
isotherms  of  interest,  Tabata  selected  isotherms  which  were 

within  the  permanent  thermocline,  assuming  that  the  thermocline 

i 
would  oscillate  as  a  unit,  although  phase  differences  between 

thermoclines  would  not  be  indicated  in  the  power  spectrum. 

Seiwell,  in  the  North  Atlantic,  chose  the  20  C  isotherm  as  being 

close  to  defining  the  top  of  the  thermocline,  or  depth  of  the 

mixed  layer. 

Defant  |l932|  in  examining  the  data  collected  by  METEOR 

anchor  stations,  used  a  similar  method,  although  he  examined 

temperature  fluctuations  at  fixed  depths.   In  order  to  compare 

his  results  with  that  of  other  investigators,  one  must  select 
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a  depth  near  the  top  of  the  thermoeline,  since  other  investigators 
used  isotherms  near  the  top  of  the  thermoeline  for  analysis. 
Decreases  in  temperature  would  indicate  the  approach  of  a  crest, 
and  increases  in  temperature  the  approach  of  a  trough. 


LaFond  and  Rao 


1954 


obtain  their  raw  data  from  the  mixed 
layer  depth,  disregarding  any  changes  above  or  below  this  depth, 
which  corresponds  to  the  top  of  the  permanent  thermoeline.   This 
would  appear  to  be  a  more  realistic  approach,  since  the  top  of 
the  thermoeline  will  approximate  a  density  discontinuity,  and 
internal  waves  should  have  their  greatest  amplitudes  at  this 
point,  and  thus  be  most  easily  detected. 

In  this  thesis,  the  fluctuations  in  depth  of  the  top,  center, 
and  bottom  of  the  thermoeline  were  used  to  investigate  internal 
waves.   The  depths  of  these  features  were  obtained  using  a  method 
originated  by  Boston  1966;  in  which  the  thermoeline  is  made  to 
approximate  a  Gaussian  distribution  of  depths  and  related 
temperatures.   The  depths  thus  obtained  were  subjected  to  a 
Fourier  analysis,  and  the  coefficients  obtained  were  harmed 
and  smoothed  over  discrete  frequency  bands. 

For  comparison  of  results  of  isotherm  analysis,  the  depths 
of  six  selected  isotherms  were  followed  through  the  time  series 
and  subjected  to  analytical  procedures  identical  to  those  used 
in  the  spectral  analysis  of  the  thermoeline. 


III.   OBTAINING  THERMOCLINE  DEPTHS 

A.   THE  GAUSSIAN  THERMOCLINE 

A  gaussian,  or  normal,  distribution  of  temperature  (T)  as  a 
function  of  elevation  (Z)  is  given  by  the  equation 


where  Z  is  positive  upwards.   The  frequency  function  corresponding 
to  this  distribution  function  is 

dTZ)  _   1  -4 


tiZ   "V2ir 


Using  a  basic  statistical  approach,  defining  a  normalized  variable 

y,  _  Z-m 


a 

where  m  is  the  mean  of  the  Z  values,  and  CT  is  the  standard 
deviation  of  the  array  of  Z  values,  we  obtain  a  new  frequency 
function  which  is 


(Z-m)2 


clT(Z)_   1  p? 
dZ   oV2ir 


This  function  has  the  following  characteristics: 

1.  It  is  symmetric  about  the  point  Z  =  m 

2.  It  has  two  symmetric  points  of  inflection  at  Z  =  m-CT 

3.  It  has  a  maximum  rate  of  change  at  Z  =  m-G"V-3 

A  change  in  the  numerical  value  of  m  causes  a  displacement  of  the 
curve  in  the  vertical  direction,  but  does  not  alter  its  form. 


10 


A  change  in  J  has  the  effect  of  altering  the  scale  in  the 
vertical  direction.   The  smaller  the  value  of   J"  ,  the  more 
concentrated  the  curve  is  about  the  .point  Z  =  m. 

If  the  water  surface  is  defined  to  be  Z  =  0,  and  the  depth  of 
the  water  considered  to  be  infinite,  then  the  function  defining 
the  vertical  temperature  distribution  is 


T(Z')  = 


1 

-( 

-oo 

D- 

(Z-mf 

n  ■ 

^J2it 

The  center  of  the  curve  T(ZT)  versus  Z  is  found  at  the  point 
Z  =  -m. 

If  this  curve  is  considered  as  a  bathythermograph  trace,  the 
following  definitions  may  be  made.   The  center  of  the  distribution 
of  T(Zf)  versus  Z  is  the  center  of  the  thermocline,  called  Z  . 
The  top  and  bottom  of  the  thermocline  will  correspond  to  the 
points  of  maximum  rate  of  change  of  slope  of  this  distribution. 
If  Z.  is  used  to  designate  the  depth  of  the  top  of  the  thermocline 
and  Z]-^  to  designate  the  depth  of  the  bottom  of  the  thermocline, 
these  definitions  become 

Zc  =  -m 

Zt  =  -m+O" 

zbt=  -m-CT 
These  features  are  illustrated  in  Figure  1A. 

To  obtain  the  necessary  values  to  apply  these  definitions, 
consider  the  first  and  second  moments  of  the  (T,Z)  curve  about 
the  origin,  Z  =  0.   The  nth  moment  of  the  curve  about  the  origin 
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n         ""tt 
is,  by  definition,  the  mean  value  of  the  array  Z  ,  called  Z  , 

where  n  is  a  positive  integer.   With  a  normal  curve,  the  first 
moment  is  identical  with  the  mean,  which  corresponds  to  the 

center  of  the  thermocline. 

J(0) 


ZC-  TpO-T(-o3)jZ(T)dT 


T(-~) 


where  T(0)  is  the  temperature  at  the  surface  and  T(-OO)  is  the 
temperature  of  the  bottom  of  the  lower  isothermal  layer  which, 
so  far,  has  been  assumed  to  be  the  bottom  of  the  ocean.   The 
second  moment 

T(0) 

1 


Z-TO)-T(-oo)|Z(T)2JT 

fr°°) 

can  be  used  to  compute  O    ,  since  the  variance  of  the  distribution. 
<J   ,  is  defined  as 

--rcr-7r  -Z^ 

In  practice,  several  approximations  to  this  exact  solution 
must  be  made.   The  temperature  T(-'JO)  is  measured  at  the  bottom 
of  the  bathythermograph  trace.   The  differentials,  dT,  become 
finite  differences,   T,  and  the  integrals  are  replaced  by 
summations . 

In  changing  from  an  integral  representation  to  a  summation 
over  finite  intervals,  a  certain  amount  of  error  in  the  value 
of     is  introduced,  due  to  the  fact  that  summation  in  this 
fashion  is  not  a  continuous  process,  whereas  integration  is. 
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To  compensate  for  this  loss  of  accuracy,  a  weighting  factor  must 


be  introduced  in  all  computations  involving  C  .   Grosfils 


1968 


calls  this  factor  sk,  and  bases  the  size  of  the  factor  on  the 
size  of  OT. 

For   ST^1.0r      sk  =  1.30 

For   St^  0.5°      sk  =  1.47 
Intermediate  values  of   OT  will  have  intermediate  values  of  sk, 
but  these  values  were  chosen  as  being  those  most  commonly 
encountered.   Using  this  factor,  the  equations  defining  the  top 
and  bottom  of  the  thermocline  become 

Zt  =  -m  +  a  sk 

Zbt=  ~m  "  ^  sk 

Using  this  approach  to  define  the  thermocline,  one  can  fully 
describe  a  thermocline  knowing  the  temperature  at  the  surface, 
the  temperature  at  the  bottom  of  the  bathythermograph  cast,  and 
depths  corresponding  to  some  constant  temperature  interval. 

The  center  of  the  Gaussian  thermocline  can  now  be  defined  as 


Zc=  Z'Tn^^- )2Tn)i 


and  since 


then 


Z(Tn)i=— 1_>  Z(Tn); 
T1  "Tbt^ 


a2-  Z(Tn)f-Z2c 

Using  these  factors,  the  top  and  the  bottom  of  the  Gaussian 
thermocline  may  readily  be  determined. 
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In  the  event  the  distribution  of  depth  with  temperature  is 
not  Gaussian,  the  inaccuracies  built  into  this  objective 
definition  are  intolerable.   Whether  or  not  this  distribution 
is  Gaussian  can  be  determined  by  computing  the  third  moment  of 
the  (T,Z)  curve  about  the  origin.   This,  however,  is  unnecessary 
when  one  applies  a  non-Gaussian  approach  to  the  thermal  structure 

B.   THE  NON-GAUSSIAN  THERMOCLINE 

The  approach  that  was  used  by  Boston  and  by  Grosfils  was 
based  on  depth  data  given  at  constant  temperature  intervals.   The 
data  used  for  this  thesis  were  temperature  data  given  at 
constant  depth  intervals.   This  difference  does  not  affect 
the  treatment  of  the  data  except  for  the  determination  of  the 
temporary  center  of  the  thermocline.   The  details  of  the  method 
used  in  this  thesis  to  find  the  center  of  the  thermocline  are 
given  in  Section  IV. 

When  the  depth  of  the  center  of  the  thermocline  is  found,  by 
whatever  method  the  data  requires,  the  data  is  then  processed  as 
two  Gaussian  thermoclines — one  based  on  the  depth  distribution  of 
temperatures  above  and  including  the  depth  Z  ,  the  other  based  on 
the  depth  distribution  of  temperatures  below  and  including  the 
depth  Zc . 

In  order  to  locate  Zt,  only  the  upper  portion  of  the  (T,Z) 
curve  is  considered.   This  portion  of  the  curve  is  "reflected" 
about  Zc,  to  produce  the  curve  shown  in  Figure  IB.   From  this 
curve       is  computed  by  methods  previously  described  and, 
using  this  value,  Z-j-  may  be  calculated. 
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FIGURE   I 
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In  order  to  determine  Zjj-j-,  a  similar  procedure  is  followed, 
except  that  in  this  case,  the  lower  portion  of  the  (T,Z)  curve 
is  "reflected"  about  Zc,  to  produce  the  curve  shown  in  Figure  1C. 
Note  that  it  is  possible  for  this  resulting  curve  to  include 
negative  depths,  but  since  the  depths  are  squared  in 
computing   CT  ^ ,  this  does  not  result  in  a  narrowing  of  the 
distribution  through  an  underestimate  of  C     . 

It  can  be  seen  that  in  the  event  the  thermocline  is  Gaussian, 
treating  it  in  this  fashion  will  not  alter  the  results  obtained 
by  treating  it  as  Gaussian.   Thus,  every  (T,Z)  input  may  be  treated 
as  a  non-Gaussian  distribution,  regardless  of  its  true  form. 
This  is  the  method  which  was  used  in  processing  the  ba thy thermogram 
data  for  this  thesis. 
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IV.   DIGITAL  METHODS  USED  TO  OBTAIN  DEPTH  INPUTS 

The  bathythermograph  data  was  obtained  from  the  National 
Oceanographic  Data  Center  in  digital  format.   For  each 
bathythermograph  trace,  the  digital  data  was  contained  on  three 
or  four  computer  data  cards,  including  the  sequential  number  of 
the  observation,  the  date  and  hour  of  the  observation, 
temperature  in  degrees  Centigrade  at  each  five  meters  of  depth, 
and  such  other  bits  of  synoptic  information  as  the  location  of 
the  observation,  the  ship  making  the  observation,  and  a  number 
identifying  the  cruise  during  which  the  data  was  taken  (Figure  2) 

The  first  phase  of  data  reduction  was  to  transfer  the 
sequential  number  of  the  observation,  the  time  and  date,  and 
the  temperatures  at  the  first  twenty  depths—which  correspond  to 
temperatures  from  the  surface  to  a  depth  of  95  meters--to  a 
computer  data  card.   This  was  not  essential  for  data  reduction, 
but  was  done  for  convenience  as  it  reduced  the  size  of  the  data 
deck  by  three  quarters.   By  inspection  of  ba thy thermograms  for 
the  entire  series  of  observations  the  permanent  thermocline  was 
in  all  cases  determined  to  be  shallower  than  95  meters  in  depth. 

The  data  deck  thus  obtained  was  then  analyzed  to  obtain  the 
depths  of  the  top,  center,  and  bottom  of  the  thermocline  using 
a  method  based  on  a  digital  computer  program  developed  by 
Grosfils  (1968   .   Grosfils'  program,  however,  was  written  for 
data  consisting  of  depths  corresponding  to  equally  spaced 
temperatures,  and  since  the  NODC  format  differs  from  this,  some 
modifications  had  to  be  made  in  the  program. 
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The  computer  first  reads  the  data  from  the  data  deck,  then 
calls  subroutine  INTERP,  which  fits  straight  lines  through  the 
given  data  points,  then  determines  the  depths  of  temperatures  at 
intervals  of  0.1°C  from  the  surface  to  95  meters.   The  temperature 
is  read  by  the  computer  as  an  integer—in  tenths  of  degrees--and 
left  in  this  form  for  the  interpolation  process,  to  eliminate 
truncation  error  in  manipulating  the  temperatures. 

In  the  program  written  by  Grosfils,  the  center  of  the 
thermocline  was  defined  as  the  point  where  the  second  difference 
of  depths  was  greater  than  zero.   Due  to  the  format  of  the  data 
used  in  his  program,  however,  a  great  deal  more  information  was 
known  about  the  thermocline  than  is  given  in  the  NODC  format. 
For  this  reason,  the  following  method  was  used  to  define  the 
center  of  the  thermocline,  utilizing  the  array  of  temperatures 
before  interpolation. 

In  subroutine  INTERP,  the  first  differences  of  temperatures 
were  calculated,  from  the  surface  (corresponding  to  T, )  downward. 
The  ith  first  difference  is  defined  as 

A^i=  T±  "  T±+1 
The  first  time  the  ith  first  difference  is  less  than  -0.3°C,  T-    is 
stored  as  ITTOP  and  the  computer  then  begins  calculating  first 
differences  again  at  a  depth  of  M-5  meters.   When  the  ith  first 
difference  then  becomes  greater  than  -0.2  C,  the  computer  stores 
Ti+1  as  ITB0T'   The  temperature  of  the  center  of  the  thermocline 
is  then  defined  as  the  arithmetic  mean  of  ITTOP  and  ITBOT,  and  is 
stored  as  ITC .   The  computer  then  interpolates  depths  and 
temperatures  as  descrived  previously. 
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Once  the  interpolation  is  completed,  subroutine  PURGE  is 
called  to  insure  that  only  one  depth  is  recorded  for  each 
temperature  value,  which  would  not  be  the  case  if  any  temperature 
appeared  more  than  once  in  the  data  which  was  read  in.   PURGE 
examines  the  interpolated  temperatures,  and  in  the  event  that 
any  temperature  is  greater  than  or  equal  to  a  temperature  at  a 
shallower  depth,  discards  the  corresponding  depth.  After 
proceeding  through  the  entire  interpolated  trace  in  this  manner, 
the  trace  is  reduced  to  one  in  which  there  is  only  one  depth 
corresponding  to  any  given  temperature,  and  it  will  be  the 
shallowest  depth  at  which  this  temperature  appears.  Although 
the  program  at  this  point  has  destroyed  all  information  about 
inversions,  this  was  not  felt  to  be  an  unreasonable  step.   The 
raw  data  were  all  examined  for  inversions  prior  to  the  interpo- 
lation step,  and  in  no  case  were  inversions  present  greater  than 
+0.1  C  in  five  meters  of  depth. 

At  this  point,  subroutine  FINDZC  is  called,  and  it  examines 
the  temperature  and  depth  information  thus  computed  and  chooses 
the  depth  of  the  first  (shallowest)  temperature  which  is  less 
than  or  equal  to  ITC.  It  stores  this  depth  as  ZC,  the  depth  of 
the  center  of  the  thermocline.  This  depth,  along  with  all  the 
remaining  temperature  and  depth  information,  is  returned  to  the 
main  program. 

Using  this  information,  the  main  program  calls  subroutine 
NONGA,  which  uses  the  procedure  described  in  Section  II  to  find 
the  top  and  bottom  of  the  thermocline  through  the  use  of 
subroutine  GAUSS1  and  GAUSS2,  respectively. 
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For  comparison  with  Tabata's  results,  the  depth  fluctuations 
of  selected  isotherms  were  also  subjected  to  Fourier  analysis. 
Depth  information  for  this  analysis  was  obtained  in  a  manner 
similar  to  that  described  above,  except  that  for  depths  of 
isotherms,  subroutine  FINDZS  is  called  from  INTERP  after  PURGE 
and  the  depths  are  determined  in  a  manner  similar  to  that  used  to 
find  ZC. 

Fluctuations  of  isotherms  and  thermoclines  with  time  are 
shown  in  Figures  3  and  4,  respectively. 
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V.   COMPUTATION  OF  POWER  SPECTRA 

The  reworked  data  obtained  as  described  in  Section  IV  were 
subjected  to  Fourier  analysis  to  produce  a  power  spectrum. 
Subroutine  RHARM,  which  is  contained  in  the  IBM  Scientific 
Subroutine  Library,  was  used  for  this  analysis.   This  subroutine 


uses  an  algorithm,  developed  by  Tukey  and  Cooley 


1965 


,  to 


transform  the  depth  fluctuations  to  a  Fourier  series  of  the  form 

v1^ 


7?  =  x   ancos  ncut+b^sin  ncot 


nrOJ 
where  HCO  is  the  frequency  of  a  given  sinusoidal  component  in 

radians  per  hour,  N  is  half  the  total  number  of  data  points  plus 

one,  t  is  time  in  hours,  and 77  is  the  instantaneous  wave  elevation 

above  its  mean  position.   The  sum  of  the  squares  of  the 

2    2 
coefficients,  an  +bn  ,  is  proportional  to  the  power  contained  in 

the  nth  component,  whose  period  is  N/(n)  hours.   aQ  is  the  mean 

value  of  the  data,  and  to  eliminate  leakage  of  this  large  value 

into  neighboring  frequency  components,  all  data  was  adjusted  to 

a  mean  value  of  zero  in  subroutine  AVERGE  before  RHARM  was  called. 

bQ  is  always  equal  to  zero,  and  is  therefore  ignored  in  subsequent 

manipulations . 

Once  the  Fourier  coefficients  were  obtained,  they  were  harmed 

by  subroutine  HANN.   This  subroutine  used  the  following  recursion 

formula  to  refine  the  coefficients: 

t 
a,  =  0.5a,  -  0.25a 

a[   =  -0.25ai_1  +  O.Sai  -  0.25ai+1  i  =  2,3,...,N-1 

aN  =_0-25aN-l  +   °-5aN 
The   sine   coefficients  were   similarly  treated.      Hanning  the   data  is 


2i+ 


a  method  of  greatly  reducing  leakage  of  power  from  frequencies 
containing  large  amounts  of  power  into  neighboring  frequencies. 
According  to  Bingham,  et.  al.  11967   if  the  Fourier  transform  is 

perfect  for^J^,  ^2»  ^3* an(i  ^  *s  none  of   the  CO  ^9  the 

leakage  without  harming  is  proportional  to   U)-  ')    and  harming 
the  coefficients  reduces  the  leakage  to  a  term  proportional 
to      Cup  Co  "  . 

The  hanned  data  was  then  plotted  on  a  graph  of  power  (spectral 
density)  versus  frequency.   Although  this  presentation  gave 
indications  of  the  frequencies  of  the  major  components  of  the 
internal  waves  present,  it  was  objectionably  noisy.   To  reduce 
this  noise,  a  method  of  power  spectrum  smoothing  was  utilized  to 
produce  the  final  spectrum.   This  method,  originated  by 
Jones  11965  was  used  in  subroutine  SMOOTH.   The  method  consists 
of  first  arranging  the  squares  of  the  non-zero  a  and  b  coef- 
ficients in  order  and  renaming  them  s,  as  follows: 

2    2    2    2    2 
a0   'al   '  1   ,a2   '2   ' 


sl   's2   's3   'SM-   's5   » 

The  smoothed  spectrum  is  obtained  by  averaging  these  numbers  in 

groups  of  2m,  except  for  the  first  group,  which  contains  m,  and 

the  last  group,  which  contains  those  numbers  which  do  not  fit 

into  groups  of  (2p  +  l)m,  where  p  is  the  number  of  groups  of 

2m  included  in  the  data. 
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S  =  N/m 
o 


m 

+  2V 


s ,  +  2  7    s . 

1    /_,   \ 


Sj_  =  N/m  y   s± 


m*1 


fcl 


Sp  =  N/mp  )  s± 
(2p-1)m*1 
where  nu  is  the  number  of  pairs  remaining  in  the  last  step  and 

N  is  the  total  number  of  Fourier  coefficients  used  in  the 

smoothing.   By  this  method  of  smoothing,  SQ  will  be  centered  at 

.25/N^t,  S!  will  be  centered  at  1/N  ht,  S2  at  2/N  5t,  and  so 

forth.   Care  must  be  taken  to  choose  m  sufficiently  large  that 

significant  smoothing  takes  place  without  choosing  it  so  large 

that  it  eliminates  detail  from  the  final  spectrum. 

Subroutine  SMOOTH  was  run  with  mTs  of  2,  3,  4,  5,  10,  20, 
and  50.   Of  these,  the  best  results  were  obtained  with  m  =  M-, 
and  this  grouping  factor  was  used  in  the  final  analysis .  After 
the  data  were  smoothed,  subroutine  PERIOD  was  called  to  generate 
counters  to  be  used  for  the  final  graphic  output. 

The  data  were  broken  into  two  groups,  the  first  containing 
the  first  512  data  points  (consecutive  numbers  1  through  512)  and 
the  second  containing  the  last  512  data  points  (consecutive 
numbers  139  through  651) .  More  confidence  may  be  placed  in 
spectral  estimates  if  the  spectral  estimates  are  the  averages 
of  several  sets  of  independent  data.   The  nearest  approximation 
to  this  which  was  possible  with  this  data  was  the  treatment  which 
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was  used.   Had  a  sample  size  of  256  points  been  chosen,  totally 
independent  data  sets  could  have  been  used,  and  confidence  in 
the  results  would  have  been  increased.   The  reduction  in  the 
number  of  data  points  by  a  factor  of  two,  however,  would  have  so 
significantly  reduced  the  accuracy  of  the  harmonic  analysis  that 
the  increased  confidence  in  the  results  would  have  been 
meaningless.   For  this  reason,  the  use  of  two  partially  independ- 
ent data  sets  of  512  points  each  was  decided  upon. 

The  two  groups  of  hanned,  smoothed  data  were  then  added 
together,  term  for  term,  by  subroutine  GROUP,  and  averaged. 
Plots  of  the  results  of  this  averaging  for  the  top,  center,  and 
bottom  of  the  thermocline  are  shown  in  Figures  5,  6,  and  7. 

Analysis  of  the  depth  fluctuations  of  selected  isotherms  was 
carried  on  in  the  same  manner  as  analysis  of  the  thermocline 
data.   The  smoothed  spectra  of  fluctuations  of  the  isotherms 
analyzed  are  shown  in  Figures  8,  9,  10,  11,  12,  and  13.   For 
comparison,  Tabata's  results  for  the  6  C  and  10  C  isotherms  are 
shown  in  Figure  1M-. 

The  complete  computer  program  used  in  data  reduction  is 
given  in  the  Computer  Program  section  of  this  thesis . 
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FIGURE     7 
Z^   SPECTRUM 
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FIGURE     8 
11°C    ISOTHERM  SPECTRUM 
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FIGURE      9 

10°C     ISOTHERM  SPECTRUM 
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FIGURE     IO 

9°C     ISOTHERM  SPECTRUM 
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FIGURE     II 
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FIGURE     12 

7°C    ISOTHERM    SPECTRUM 
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FIGURE     I  3 

6°C    ISOTHERM  SPECTRUM 


36 


LN  FREQUENCY(hr'x128) 
FIGURE     14 

6°C&10°C  ISOTHERM  SPECTRA  (tabata,1965) 
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VI.   RESULTS 

In  his  paper,  Tabata  reports  a  strong  semidiurnal  component 
in  the  internal  waves  analyzed  through  depth  fluctuations  of  the 
6CC  through  10  C  isotherms,  with  greater  power  in  higher 
frequencies  noted  in  analysis  of  the  deeper  isotherms.   This 
general  behavior  was  also  noted  in  the  digitally  analyzed  data 
used  in  this  thesis . 

A  power  peak  centered  at  an  11.64  hour  period  was  present 
in  the  spectra  for  all  the  isotherms  analyzed,  as  well  as  in  the 
spectra  of  the  top,  center,  and  bottom  of  the  thermocline — the 
shallowest  feature  studied—also  shows  a  peak  corresponding  to  a 
period  of  25.6  hours,  probably  associated  with  a  principal  lunar 
diurnal  (OjJ)  tidal  component. 

The  power  spectrum  due  to  fluctuations  of  the  bottom  of  the 
thermocline  indicates  a  shift  of  power  to  higher  frequencies, 
with  peaks  associated  with  periods  of  approximately  4.5  and  3.5 
hours,  and  a  maximum  energy  associated  with  waves  of  9.14  hour 
period.   Tabata  reports  waves  of  periods  equal  to  five  hours 
and  slightly  less  than  four  hours  present  for  the  deepest  of 
the  isotherms  which  he  studied,  with  which  these  data  agree.  In 
no  case,  however,  does  he  report  waves  of  periods  near  nine 
hours.   Relative  peaks  in  spectral  density  having  periods  between 
nine  and  ten  hours  (most  often  at  a  period  of  9.85  hours)  were 
found  associated  with  all  features.   Waves  with  periods  of 
16.00  hours,  which  is  very  close  to  the  inertial ': period  lof  15.7 
hours  at  Ocean  Station  "P"  were  also  noted  in  all  spectra 
produced. 
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In  the  original  plan  for  the  analysis,  the  fluctuations  of  the 
thickness  of  the  thermocline  were  also  to  be  investigated,  to 
determine  whether  the  thermocline  was  oscillating  in  the  funda- 
mental mode — in  which  the  entire  thermocline  oscillates  as  a 
unit — or  mode  two — in  which  the  top  and  bottom  of  the  thermocline 
are  180  out  of  phase,  resulting  in  a  periodic  thickening  and 
thinning  of  the  thermocline fs  vertical  extent.   The  spectra  pro- 
duced by  analysis  of  the  thickness  showed  nearly  constant  spec- 
tral density  throughout  the  entire  range  of  frequencies  studied, 
except  for  large  contributions  in  the  low  frequency  end  of  the 
spectrum.   This  suggested  that  the  thermocline  was  oscillating  in 
the  fundamental  mode,  and  the  spectral  densities  at  high 
frequencies  were  due  to  noise.   To  verify  this,  the  correlations 
of  the  top,  center,  and  bottom  of  the  thermocline  were  computed. 
It  was  found  that  the  depth  fluctuations  of  each  of  these 
features  with  the  others  was  greater  than  +.90  which  confirms 
the  fact  that  the  entire  thermocline  was  oscillating  as  a  unit. 
This  conclusion  agrees  with  the  results  published  by  Tabata,  who 
found  phase  differences  on  the  order  of  -  20°  between  oscillations 
of  the  isotherms  he  analyzed. 

The  spectrum  associated  with  the  top  of  the  thermocline 
showed  a  very  large  contribution  due  to  waves  of  very  long 
period.   Examination  of  Figure  3  shows  the  reason  for  this.   The 
top  of  the  thermocline  was  slowly  sinking  with  time,  presumably 
due  to  mixing  in  the  upper  layer.   This  effect,  although  small 
in  magnitude,  would  have  been  identified  in  the  analysis  as  an 
internal  wave  with  a  period  equal  to  the  total  length  of  the 
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record,  which  in  this  case  was  512  hours.   In  addition,  during 
the  period  16  August  to  21  August,  the  water  column,  which  was 
almost  perfectly  isothermal  above  and  below  the  thermocline  at 
all  other  times,  developed  near-surface  transient  the rmoc line s . 
As  the  non-Gaussian  analysis  used  to  obtain  thermocline  depths 
did  not  eliminate  transient  thermoclines  before  computing  the 
depths  associated  with  the  main  thermocline,  these  near-surface 
transients  resulted  in  the  calculation  of  an  abnormally  large 
variance  by  subroutine  GAUSS1,  which  had  the  effect  of  moving 
the  calculated  top  of  the  thermocline  closer  to  the  surface  than 
it  actually  was.   The  fact  that  the  time  period  during  which  this 
occurred  was  included  in  both  sets  of  512  hourly  observations 
meant  that  an  additional  large  power  component  associated  with 
a  period  of  512  hours  was  found  in  the  analysis.   This  is  an 
apparent  effect,  rather  than  a  real  one,  and  even  when  combined 
with  the  long-term  sinking  of  the  thermocline,  fails  to  obscure 
the  results  in  the  frequency  range  of  interest. 

Examination  of  Figure  3  also  reveals  that  the  center  of  the 
thermocline  also  showed  a  slight  downward  displacement  with  time, 
but  of  less  magnitude  than  that  of  the  top  of  the  thermocline. 
This  resulted  in  a  relatively  large  amount  of  power  associated 
with  long-period  internal  waves  at  the  center  of  the  thermo- 
cline, but  much  smaller  than  the  power  contribution  to  the 
spectrum  of  the  top  of  the  thermocline.   The  bottom  of  the 
thermocline  showed  no  change  in  mean  depth  with  time,  which 
indicates  that  surface  layer  mixing  extended  only  to  a  depth 
slightly  below  the  center  of  the  thermocline .   Although  in  most 
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cases  large  amounts  of  power  were  associated  with  low  frequencies, 
these  frequencies  were  not  considered  of  interest,  and  are  not 
indicated  on  the  graphs  of  the  final  spectra.   A  complete  presen- 
tation of  all  data  produced  by  the  analysis  is  presented  in  the 
COMPUTER  OUTPUT  section. 

Tabata's  results  for  the  6°C  and  10°  C  isotherms  are  presented 
in  Figure  14  in  a  semilog  plot  for  comparison  with  the  results  of 
the  digital  analysis.   Although  these  curves  differ  in  detail 
from  those  presented  for  the  same  isotherms  analyzed  digitally  in 
this  thesis,  the  overall  form  is  strikingly  similar.   For  both 
analyses  of  the  10°C  isotherm,  the  greatest  peak  spectral  density 
is  associated  with  waves  whose  period  is  greater  than  the  semi- 
diurnal tidal  period.   For  both  analyses  of  the  6°C  isotherm, 
the  spectral  density  is  rather  well-distributed  throughout  the 
higher  frequency  components,  which  also  have  a  greater  mean  spec- 
tral density  than  the  higher  frequency  components  of  the  analyzed 
10  C  isotherm  data.   This  indicates  that  the  higher  frequency 
components  showed  an  increase  in  amplitude  with  increasing  depth, 
a  behavior  which  was  also  noted  by  Haurwitz,  Stommel,  and  Munk 
1959]  . 

The  Brunt-Vaisala  frequency  associated  with  the  thermal 

2      1 

structure  of  this  time  series  is  3.79  x  1Q~^  sec  .  This  indi- 
cates that  the  thermal  structure  is  capable  of  supporting  inter- 
nal waves  of  any  period  greater  than  26.4  seconds.  The  Nyquist 
frequency  for  a  digital  time  series  sampled  at  hourly  intervals 
is  0.5  hr  .  Thus,  there  is  a  large  range  of  frequencies  which 
may  be  present  but  of  which  the  data  does  not  permit  investigation. 
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The  causes  of  the  internal  waves  with  periods  near  2M-,  16, 
or  12  hours  seem  obvious .   The  24-  and  12  hour  periods  are 
probably  associated  with  internal  waves  which  are  being  driven 
by  diurnal  or  semidiurnal  tidal  forces.  As  mentioned  previously, 
the  16  hour  period  seems  to  be  associated  with  f ne rtial  iintarnal .'. 
waves.   The  causes  for  the  internal  waves  with  periods  between 
four  and  six  hours — indicated  both  in  this  thesis  and  in  TabataTs 
study — are  not  so  easily  theorized. 

Another  possible  source  of  internal  waves  is  baroclinic 
Rossby  waves,  but  this  type  of  wave  is  associated  with  extremely 
low  frequencies,  with  periods  on  the  order  of  days,  so  it  would 
seem  that  these  waves  will  not  explain  the  higher  frequency 
waves  indicated  by  the  analysis .   Two  other  possible  causes  for 
these  waves  do  exist,  however.   They  might  be  long-period  waves 
which  have  been  reflected  from  a  coast  and  have  interacted  with 
varying  density  gradients,  which  act  as  a  weak  filter,  but  this 
theory  cannot  be  tested  unless  one  knows  the  direction  of 
propagation  of  the  higher  frequency  internal  waves . 

The  other  possible  source  of  these  waves  is  a  non-linear 
interaction  between  surface  pressure .disturbances,  but  in  order 
to  examine  this  as  a  possible  cause,  one  must  have  records  of 
both  surface  waves  and  weather.  Although  weather  records  exist 
for  the  time  during  which  the  bathythermographs  were  taken,  no 
surface  wave  information  is  available.   To  be  of  value  in 
investigating  this  process  as  a  causative  agent,  the  surface 
wave  information  would  have  to  be  very  detailed,  and  wave  records 
of  the  required  detail  can  only  be  obtained  in  the  deep  ocean 
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through  highly  sophisticated  and  careful  measurements  which  are 
not  compatible  with  the  mission  of  a  weather  ship. 

In  conclusion,  it  seems  that  little  can  be  said  about  these 
higher-frequency  internal  waves  except  that  they  exist.   In  order 
to  say  more,  one  would  require  a  great  deal  more  data  than  is 
available.   Little  has  been  done  in  the  investigation  of  internal 
waves  to  date,  and  one  can  see  that  it  represents  a  fertile 
subject  for  further  research. 
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VII.   RECOMMENDATIONS 

During  the  execution  of  this  thesis,  a  number  of  changes  in 
the  format  and  amount  of  data  were  wished.   Had  these  changes 
been  made,  the  analysis  of  the  data  would  have  been  simplified 
and  the  accuracy  of  the  results  enhanced.   For  these  reasons,  the 
following  recommendations  are  made . 

In  digitizing  mechanical  bathythermograph  traces,  it  would 
greatly  assist  the  user  of  the  digitized  data  if  the  data 
through  a  thermocline  were  recorded  as  depths  at  equally  spaced 
temperatures,  instead  of  temperatures  at  equally  spaced  depths 
as  is  now  done.   This  would  provide  an  increased  number  of  points 
in  the  region  of  maximum  gradient,  enabling  the  digitized  data 
to  be  used  to  more  accurately  reconstruct  the  original  trace. 
The  method  presently  used  to  digitize  expendable  bathythermograph 
traces,  which  records  a  data  point  at  the  location  where  the 
trace  departs  from  linearity  by  a  specified  amount,  allows  use 
of  digitized  data  to  almost  perfectly  reconstruct  the  trace. 

If  salinity  data  were  also  available  for  analysis,  it  would 
be  a  simple  task  to  compute  depths  of  isopycnals  using  a 
computer.   If  terms  to  second  order  involving  salinity  and 
pressure  were  introduced  to  the  already-known  temperature 
information,  depth-density  profiles  of  much  greater  accuracy 
would  result,  and  this  in  turn  could  be  analyzed  to  produce 
internal  wave  information  of  higher  quality. 

Closer  spacing  of  observations  would  enable  analysis  of 
higher  frequency  components,  since  through  aliasing  the  smallest 
wave  period  that  can  be  detected  using  a  sampling  interval  of  one 
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hour  is  two  hours.   Also,  as  was  suggested  by  Tabata,  utiliz i  i 
data  presently  available  over  extremely  long  periods,  with  t  as 
Large  as  a  year,  would  allow  analysis  of  extremely  low  frequency 
internal  waves. 

As  was  indicated  in  the  Results  section  of  this  thesis, 
there  is  a  need  for  information  regarding  the  direction  of 
propagation  of  the  internal  waves  indicated.   This  could  be  got- 
ten using  a  cross  array,  consisting  of  two  orthogonal  arrays 
oi  the  type  described  by  Gilchrist  il966 j  for  use  in  investi- 
gating the  direction  of  propagation  of  surface  waves.   SeveraJ 
such  arrays  would  be  required  to  allow  investigation  of  a  large 
range  of  possible  wavelengths,  but  it  would  seem  that  the 
information  obtained  as  a  result  of  such  an  expenditure  of  effort 
and  equipment  would  be  a  highly  significant  contribution  to  the 
present  knowledge  of  internal  waves. 
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APPENDIX  A 
EXAMPLE  OF  COMPUTATION 

The  following  four  figures  indicate  graphically  the  operation 
of  the  subroutines  used  in  computing  the  depths  used  in  harmonica 
analysis.   The  data  is  from  bathythermograph  number  392  (Fig.  2).- 

Figure  15  indicates  the  operation  of  subroutine  INTERP  in  fit- 
ting straight  line  segments  between  raw  data  points,  and  the 
choosing  of  a  temperature  corresponding  to  the  center  of  ther- 
mocline  (ITC) . 

Figure  16  shows  the  operation  of  subroutine  PURGE  in  elimi- 
nating small  inversions  and  duplications  of  temperatures,  reducing 
the  data  to  a  set  of  temperature -depth  pairs  from  the  surface  to 
the  bottom  of  the  trace  at  an  interval  of  0.1°C.   This  figure 
also  indicates  the  operation  of  subroutines  FINDZC  and  FINDZS  in 
determining  the  depth  of  the  center  of  the  non-Gaussian  thermo- 
cline  and  the  depths  of  the  selected  isotherms. 

Figures  17  and  18  show  the  thermal  structures  used  in  GAUSS1 
and  GAUSS2  to  obtain  the  depth  of  the  top  and  bottom  of  the 
thermocline,  respectively. 
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FIGURE  15 
Subroutine    INTERP 
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FIGURE  16 
Subroutine  PURGE, 

Subroutine  FINDZQ 

and 
Subroutine  FINDZS 


48 


5 


7 


TE  MPERATURE  (  °C) 
8  9         10       11         12 

T 


I 
I— 
Q_ 
LLl 
Q 


FIGURE  17 
Subroutine  GAUSS1 
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FIGURE   18 
Subroutine  GAUSS2 
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9.54 

39.34 

40.22 

3.37 

2.46 

8.65 

8.75 

3.28 

8.32 

6.48 

8.51 

3.20 

7.52 

13.95 

11.41 

3.12 

10.27 

7.18 

16.97 

3.05 

18.  60 

16.51 

13.  C8 

2.98 

9. 50 

9.84 

18.43 

2.91 

5.25 

5.91 

5.05 

2.84 

5.94 

10.51 

9.47 

2.78 

8.10 

15.46 

18.04 

2.72 

6.22 

9.52 

15. C5 

2.67 

12.63 

9.96 

9.56 

2.61 

2.99 

9.36 

4.56 

2.56 

7.37 

6.54 

8.84 

2.51 

13.50 

12.56 

36.85 

2.46 

17.12 

13.96 

12.55 

2.42 

6.24 

12.82 

22.43 

2.37 

4.75 

13.57 

7.89 

2.33 

1.19 

20.42 

35.32 

2.29 

7.4C 

14.15 

7.85 

2.25 

4.74 

9.45 

4.69 

2.21 

8.02 

8.11 

12.78 

2.17 

9.94 

7.15 

1.35 

2.13 

15.54 

7.77 

11.75 

2.10 

18.48 

14.44 

11.56 

2.06 

0.76 

4.66 

13.29 

2.03 

0.82 

4.41 

5.74 

1.99 

3.93 

3.27 

28.37 
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SPECTRAL 

ESTIMATES 

FOR    GROUP 

NUMBER    2 

PERIOD 

11    DEG 

10    DEG 

9    DEG 

2b6.CC 

52.16 

69.90 

108.77 

128. CC 

26. CI 

31.23 

38.58 

64. CC 

54. 4C 

66.58 

75.33 

42.67 

27.93 

32.99 

36.67 

32. CC 

12.73 

13.47 

14.  18 

25. 6C 

47.71 

44.67 

42.49 

21.32 

21. 5C 

21.45 

25.32 

18.29 

35.28 

32.62 

24.70 

16. CC 

40.05 

45.20 

54.94 

14.22 

8.69 

13.78 

18.46 

12.  80 

13.88 

12.51 

8.94 

1 1.64 

60.50 

53.13 

52.29 

10.67 

31.83 

29.07 

3C.63 

9.85 

74.73 

75.01 

68. C6 

9.14 

62.45 

62.82 

55.57 

8.53 

46.58 

43.35 

40.45 

8.00 

36.11 

33.86 

33.10 

7.53 

35.29 

35.05 

32.63 

7.11 

31  .00 

26.05 

24. tO 

6.74 

18.77 

20.93 

22.63 

6.4C 

29.29 

28.22 

22.97 

6.1C 

3.03 

4.18 

3.45 

5.82 

57.95 

55.82 

51.14 

5.57 

10.39 

12.71 

14.96 

5.33 

13.85 

16.17 

18.22 

5.12 

28. ie 

27.08 

34.45 

4.92 

28.18 

28.90 

28.48 

4.74 

45.51 

49.04 

41  .54 

4.57 

9.48 

11.09 

13.60 

4.41 

24. oe 

25.39 

24.76 

4.27 

12.55 

14.16 

15.29 

4.  13 

10.37 

11.26 

1C.C1 

4.00 

22.87 

21.64 

17.14 

3.88 

15.23 

13.69 

11.97 

3.76 

18.28 

19.52 

2C.96 

3.66 

8.13 

7.64 

5.53 

3.56 

9.43 

10.50 

13.65 

3.46 

19.05 

24.26 

32.10 

3.37 

9.19 

8.84 

9.C3 

3.28 

8.85 

7.8C 

8.40 

3.20 

14.66 

15.48 

14.16 

3.12 

8.10 

7.28 

5.28 

3.C5 

14.92 

17.01 

16.74 

2.98 

9.45 

7.67 

8.27 

2.91 

11.37 

9.91 

11.06 

2.84 

6.44 

7.14 

1C.95 

2.78 

13.31 

13.44 

12.51 

2.72 

9.73 

9.70 

9.37 

2.67 

9.18 

8.74 

9.78 

2.61 

12.80 

14.21 

16.39 

2.56 

2.42 

4.57 

8.24 

2.51 

8.99 

9.65 

14.19 

2.46 

8.45 

9.08 

11.27 

2.42 

11.11 

12.28 

13.77 

2.37 

5.68 

8.81 

9.51 

2.33 

13.06 

15.17 

16.70 

2.29 

9.79 

10.11 

11.41 

2.25 

3.36 

3.40 

7.66 

2.21 

13.41 

13.17 

13. CO 

2.17 

2.21 

2.81 

6.88 

2.13 

4.79 

4.42 

5.49 

2.  1C 

12.13 

12.21 

9.89 

2.C6 

7.11 

8.91 

9.91 

2.03 

4.67 

4.57 

4.76 

1.9  9 

4.76 

4.45 

3.54 
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SPECTRAL 

ESTIMATES 

FOR    GROUP 

NUMBER    2 

PERIOD 

8    DEG 

7    DEG 

6    DEG 

256.00 

194.69 

213.73 

173.83 

128.00 

57.58 

129.74 

147.77 

64.00 

76.76 

133.56 

193.10 

42.67 

55.49 

100.73 

79.06 

32. OC 

21.23 

76. 6C 

65.44 

25.60 

42.05 

55.35 

18.11 

21.33 

21.95 

36.54 

40.63 

18.29 

11.65 

30.71 

25.92 

16. CO 

46.29 

66.71 

135.36 

14.22 

18.27 

19. 70 

42.29 

12.80 

4.94 

6.52 

50.43 

11.64 

48. CO 

53.58 

48.18 

10.67 

38.98 

49.28 

38.86 

9.85 

63.58 

68.84 

55.95 

9.14 

55.29 

80.36 

57.35 

8.53 

39.34 

31.65 

7.71 

8.00 

33.78 

43.01 

6C.06 

7.53 

30.41 

39.84 

28.44 

7.11 

17.66 

54.21 

116.32 

6.74 

23.74 

23.54 

15.47 

6.40 

20.10 

21.76 

36.51 

6.10 

9.99 

6.70 

4.73 

5.62 

71.99 

61.25 

52.81 

5.57 

15.34 

36.78 

23.22 

5.33 

9.67 

6.59 

17.98 

5.12 

34.92 

35.62 

21.27 

4.92 

30.63 

28.13 

28.39 

4.74 

39.80 

23.08 

32.17 

4.57 

12.48 

47.87 

64.83 

4.41 

21.42 

65.73 

38.23 

4.27 

15.05 

26.20 

10.77 

4.13 

6.34 

18.78 

14.18 

4.00 

15.61 

14.13 

20.43 

3.88 

12.49 

15.41 

23.33 

3.76 

14.23 

30.86 

51.93 

3.66 

6.18 

16.89 

18.91 

3.56 

13.21 

33.64 

70.53 

3.46 

33.28 

42.80 

33.14 

3.37 

11.39 

11.50 

49.13 

3.28 

14.58 

15.97 

19. C8 

3.20 

9.39 

21.56 

20.77 

3.12 

9.07 

25.88 

15.84 

3.05 

15.49 

22.13 

19.59 

2.98 

10.36 

6.02 

9.61 

2.91 

11.80 

28.53 

43.65 

2.84 

7.61 

15.23 

18.90 

2.78 

8.62 

10.36 

20.08 

2.72 

11.20 

29.23 

23.79 

2.67 

18.19 

16.42 

13.42 

2.61 

11.65 

32.13 

48. CI 

2.56 

7.39 

10.89 

22.35 

2.51 

17.48 

22.07 

56.28 

2.46 

10.91 

28.46 

45.81 

2.42 

21.95 

48.08 

41.00 

2.37 

8.38 

7.11 

21.15 

2.33 

18.11 

35.27 

19.77 

2.29 

12.87 

10.35 

17.23 

2.25 

17.15 

9.62 

11.68 

2.21 

19.65 

12.33 

22.57 

2.17 

7.90 

18.77 

15.58 

2.13 

7.50 

10.18 

35.26 

2.10 

8.71 

4.04 

12.44 

2.06 

9.69 

14.64 

8.94 

2.03 

5.92 

9.51 

3.79 

1.99 

3.34 

19.49 

14.29 
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AVERAGEC 

SPECTRAL 

ESTIMATES 

NR 

PERICC 

TCP 

CENTER 

BCTTCV 

0.5 

256. OC 

1611.60 

113.37 

24. 2C 

L.O 

128. CO 

1166.87 

55.85 

54.93 

2.G 

64.00 

91.42 

83.71 

27.22 

3.0 

42.67 

64.69 

29.58 

19.80 

4.0 

32. CC 

47.53 

20.40 

14.77 

5.C 

25.60 

110.29 

50.83 

23.66 

6.C 

21.33 

80.94 

28.71 

48.73 

7.C 

18.29 

44.66 

33.89 

32.26 

8.C 

16. CC 

40.23 

64.42 

70.67 

9.0 

14.22 

26.10 

30.40 

22.52 

10.0 

12.80 

30.06 

26.59 

33.76 

11.0 

11.64 

96.46 

109.93 

68.12 

12.0 

1C.67 

33.56 

46.61 

28.24 

13.0 

9.85 

43.40 

93.51 

36.28 

14. C 

9.  14 

27.96 

55.63 

87.34 

15. C 

8.53 

16.37 

57.95 

26.39 

16.0 

8.00 

13.02 

29.93 

39.37 

17. C 

7.53 

29.98 

43.93 

30.50 

18. C 

7.11 

32.42 

32.37 

45.84 

19.0 

6.74 

12.29 

33.62 

11.73 

2C.C 

6.4C 

47.12 

25.09 

5.10 

21.0 

6.10 

11.17 

6.83 

4.01 

22.0 

5.82 

11.60 

55.35 

23.84 

23. C 

5.57 

7.52 

18.66 

12.02 

24.0 

5.33 

19.73 

31.02 

17.48 

25. C 

5.  12 

11.21 

23.54 

10.10 

26. C 

4.92 

14.69 

22.72 

14.94 

27. C 

4.74 

25.99 

44.32 

24.21 

28.0 

4.57 

22.62 

24.08 

43.84 

29.0 

4.41 

22.06 

42.40 

26.43 

30.0 

4.27 

6.30 

15.57 

10.43 

31.  C 

4.  13 

7.85 

14.23 

20.37 

32.0 

4.00 

12.71 

23.95 

23.87 

33.0 

3.86 

11.18 

15.63 

29.13 

34.0 

3.76 

10.30 

14.22 

17.35 

35.0 

3.66 

5.86 

6.54 

24.05 

36. C 

3.56 

4.42 

12.81 

24.06 

37.0 

3.46 

10.33 

29.11 

43.42 

38. C 

3.37 

3.46 

13.47 

8.18 

39.0 

3.26 

6.93 

8.48 

8.11 

40. 0 

3.20 

11.78 

19.17 

16.63 

41.0 

3.12 

6.28 

7.76 

24.15 

42. C 

3. 05 

23.56 

11.92 

19.51 

43.0 

2.98 

15.56 

13.86 

1C.94 

44.0 

2.91 

6.93 

13.24 

8.39 

45.0 

2.84 

11.75 

8.46 

12.39 

46. C 

2.78 

7.45 

17.10 

22. 3C 

47. C 

2.72 

10.31 

12.39 

9.96 

48. C 

2.67 

10.40 

12.31 

8.61 

49.0 

2.61 

5.14 

10.94 

3.61 

5C.C 

2.56 

9.26 

12.73 

7.58 

51.0 

2.51 

11.91 

10.45 

22.76 

52. C 

2.46 

15.58 

10.50 

12.08 

53.  C 

2.42 

6.92 

13.79 

18.09 

54.0 

2.37 

3.95 

7.77 

9.58 

55.0 

2.33 

1.79 

13.99 

21.94 

56. C 

2.29 

6.64 

15.37 

14.00 

57.0 

2.25 

9.22 

12.06 

6.94 

58. C 

2.21 

7.71 

10.69 

10.59 

59. C 

2.17 

9.89 

9.37 

6.8C 

60. C 

2.13 

15.85 

8.97 

10.22 

61. C 

2.1C 

14.94 

15.86 

15.74 

62.0 

2.06 

1.21 

5.88 

12.63 

63.0 

2.03 

1.48 

5.70 

9.28 

64.3 

1.99 

3.29 

2.45 

26.72 
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AVERAGEC 

SPECTRAL 

ESTIMATES 

NR 

PERICD 

11    DEG 

1C    OEG 

9    DEG 

0.5 

256. OC 

162.91 

165.76 

159.68 

1.0 

128. OC 

44.9  2 

51.19 

65.29 

2.0 

64.00 

67.77 

82.38 

94.57 

3.0 

42.67 

21.71 

25.75 

29.99 

4.0 

32.00 

19.92 

21.37 

22.67 

5.0 

25.60 

56.58 

54.65 

52.24 

6.C 

21.33 

25.39 

23.70 

23.31 

7.C 

18.29 

40.2  7 

37.08 

31.31 

8.0 

16.00 

69.02 

69.30 

76.93 

9.C 

14.22 

31.83 

34.73 

38.67 

10.0 

12.80 

25.37 

26.21 

25.89 

11. C 

11.64 

115.22 

107.12 

104.43 

12.0 

10.67 

45.80 

43.77 

45.55 

13.0 

9.85 

80.59 

82.38 

8C.97 

14. C 

9.14 

47.02 

47.02 

42.24 

15.0 

8.53 

56.39 

54.23 

53.42 

16. C 

8.00 

31.55 

28.40 

27.25 

17.0 

7.53 

39.24 

38.39 

36.31 

18.0 

7.11 

33.06 

28.83 

26.87 

19.0 

6.74 

35.61 

37.90 

40.03 

20. C 

6.4C 

24.61 

24.21 

20.84 

21.0 

6.10 

6.70 

7.14 

5.89 

22.0 

5.82 

53.21 

51.83 

46.86 

23.0 

5.57 

15.05 

16.19 

16.22 

24.0 

5.33 

17.53 

20.50 

23. 7C 

25. C 

5.12 

20.67 

20.47 

25.93 

26.0 

4.92 

26.93 

27.85 

28.97 

27. C 

4.74 

42.42 

43.27 

37.00 

28. C 

4.57 

12.75 

14.84 

16.32 

29.0 

4.41 

35.58 

38.27 

38.00 

30.0 

4.27 

12.98 

13.80 

14.15 

31.0 

4.13 

17.97 

17.89 

16.36 

32.0 

4.00 

22.91 

22.90 

20.63 

33.0 

3.88 

15.59 

14.99 

13.49 

34.0 

3.76 

13.94 

14.86 

15.79 

35.0 

3.66 

9.60 

9.15 

7.98 

36. C 

3.56 

13.33 

14.41 

17.13 

37.0 

3.46 

18.46 

22.28 

25.80 

38.0 

3.37 

16.06 

16.02 

13.46 

39. C 

3.28 

10.43 

9.69 

8.32 

40.0 

3.20 

15.16 

16.00 

14.42 

41.0 

3.12 

6.59 

6.16 

4.83 

42.0 

3.05 

10.39 

10.86 

11.14 

43.0 

2.98 

13.66 

13.17 

14.02 

44.0 

2.91 

15.85 

14.82 

14.36 

45.0 

2.84 

8.07 

7.76 

9.75 

46.0 

2.78 

14.55 

15.12 

13.94 

47. C 

2.72 

11.78 

11.52 

12.65 

48.0 

2.67 

10.37 

10.05 

10.09 

49.0 

2.61 

13.77 

13.65 

14.83 

50.0 

2.56 

7.85 

10.26 

14.12 

51.0 

2.51 

8.62 

8.10 

10.24 

52.0 

2.46 

10.28 

9.61 

10.20 

53.0 

2.42 

13.99 

13.90 

14.68 

54.0 

2.37 

4.35 

5.72 

5.93 

55.0 

2.33 

10.87 

11.29 

11.48 

56. C 

2.29 

8.66 

9.62 

11.64 

57.0 

2.25 

6.74 

7.57 

11.46 

58.0 

2.21 

13.78 

13.64 

13.43 

59.0 

2.17 

4.88 

5.52 

9.06 

60.0 

2.13 

5.27 

5.15 

6.04 

61.0 

2.10 

11.61 

11.83 

11.01 

62.0 

2.06 

5.02 

6.65 

8.12 

63.0 

2.03 

4.61 

4.95 

4.85 

64.3 

1.99 

3.40 

3.35 

3.83 
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AVERAGEC 

SPECTRAL 

ESTIMATES 

NR 

PERIOD 

8    DEG 

7    DEG 

6    DEG 

C.5 

256. CC 

179.73 

196.59 

2C4.59 

l.C 

128. 0C 

93.43 

177.30 

2C9.18 

2.0 

64.00 

95.89 

148.00 

182.67 

3.C 

42.67 

44.80 

74.02 

91.57 

4.C 

32. CC 

32.45 

65.40 

64.49 

5.C 

25.60 

52.30 

68.61 

48.65 

6.C 

21.33 

14.06 

47.53 

44.27 

7.C 

18.29 

20.34 

68.38 

52.25 

8.C 

16. OC 

76.66 

83.28 

124.55 

9.C 

14.22 

48.05 

45.58 

53.42 

10. C 

12.80 

28.  13 

45.77 

91.  1C 

11.  C 

11.64 

99.13 

1C8.10 

1C8.25 

12. C 

10.67 

54.81 

65.62 

49.64 

13.0 

9.85 

73. 9C 

80.53 

5e.43 

14. C 

9.14 

41.56 

69.44 

53.3  8 

15. C 

8.53 

58.27 

41.26 

14.81 

16.0 

8.00 

30.92 

42.27 

67.63 

17.0 

7.53 

32.71 

43.83 

54.71 

18. C 

7.11 

20.83 

46.77 

1C1.24 

19.0 

6.74 

40.58 

33.14 

22.06 

2C.C 

6.4C 

20.8  8 

16.98 

28.16 

21.0 

6.10 

7.49 

7.14 

9.53 

22. C 

5.82 

66.69 

82.23 

73.71 

23.0 

5.57 

16.37 

36.31 

3C.40 

24.0 

5.33 

13.82 

22.67 

39.86 

25. C 

5.12 

28.61 

22.61 

13.95 

26. C 

4.92 

32.86 

34.97 

29.38 

27.0 

4.74 

35.15 

20.73 

28.57 

28. C 

4.57 

13.02 

38.68 

55.47 

29.0 

4.41 

37.33 

53.35 

43.94 

30.0 

4.27 

13.75 

22.07 

21.33 

31. C 

4.13 

16.15 

30.67 

38.65 

32.0 

4.CC 

19.39 

21.20 

3C.8C 

33.0 

3.88 

11.61 

15.19 

25.82 

34. C 

3.76 

12.60 

31.87 

48.79 

35.0 

3.66 

11.67 

15.52 

3C.17 

36. C 

3.56 

16.64 

32.46 

82.58 

37.0 

3.46 

37.94 

64.57 

67.34 

38.0 

3.37 

16.01 

10.82 

31.14 

39. C 

3.28 

10.31 

19.14 

13.02 

40. C 

3.20 

10.90 

18.52 

27.46 

41.0 

3.12 

6.63 

27.26 

21.08 

42. C 

3.05 

12.27 

22.53 

23.33 

43.0 

2.98 

14.78 

8.36 

1C.4C 

44.  C 

2.91 

16.91 

27.62 

31.91 

45. C 

2.84 

9.46 

19.44 

43.64 

46. C 

2.78 

10.85 

14.34 

16.11 

47. C 

2.72 

16.09 

27.71 

24.32 

48. C 

2.67 

13.39 

17.59 

15. 8C 

49. C 

2.61 

13.50 

23.79 

28.43 

50. C 

2.56 

12.83 

14.74 

15.59 

51.0 

2.51 

14.42 

21.48 

39.87 

52.0 

2.46 

7.06 

23.96 

35.20 

53.0 

2.42 

20.16 

39.66 

44.06 

54. C 

2.37 

6.14 

9.20 

30.64 

55. C 

2.33 

11.93 

25.52 

14.20 

56. C 

2.29 

12.88 

11.92 

16.34 

57. C 

2.25 

18.60 

18.45 

19.36 

58. C 

2.21 

16.77 

13.02 

23.25 

59. C 

2.17 

11.35 

26.15 

28.85 

60. C 

2.13 

7.22 

14.51 

42.62 

61. C 

2.1C 

13.29 

9.18 

9.39 

62.0 

2.06 

11.97 

25.83 

11.03 

63. C 

2.03 

6.42 

10.37 

3.68 

64.3 

1.99 

4.38 

24.17 

15.11 
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C  MAIN  PROGRAM 

C  THIS  PROGRAM  READS  TEMPERATURES  AT  EVEN  DEPTH  INCREMENTS 

C  OF  FIVE  METERS,  CALLS  SUBROUTINES  TO  FIND  DEPTHS  OF  THE 

C  TOP,  CENTER,  AND  BOTTOM  OF  THE  THERMOCLINE,  AND  THE 

C  DEPTH  OF  SIX  SELECTED  ISOTHERMS.   THESE  DEPTHS  ARE  SUB- 

C  JECTED  TO  A  FAST  FOURIER  ANALYSIS  TO  DETERMINE  THEIR 

C  FOURIER  COEFFICIENTS,  WHICH  ARE  THEN  HANNED,  SMOOTHED  IN 

C  FREQUENCY  BANDS,  AND  STORED  FOR  EACH  GROUP  OF  2**N  DATA 

C  CARDS.   THE  PROGRAM  MAY  ALSO  BE  USED  TO  PRODUCE  GRAPHIC 

C  DISPLAYS  OF  FLUCTUATIONS  OF  DEPTHS  WITH  TIME,  AND  THE 

C  FINAL  AVERAGED  POWER  SPECTRA. 

C 

C  VARIABLES  USED  IN  THIS  PROGRAM  ARE — 

C  I  GROUP— THE  NUMBER  OF  GROUPS  OF  2**N  DATA  CARDS  USED 

C  IN  EACH  ANALYSIS 

C 

C  IT— THE  ARRAY  OF  TEMPERATURES  IN  EACH  TRACE  AS  READ 

C  FROM  DATA  CARDS 

C 

C  Z11,Z10,Z09,Z08,Z07,Z06 — THE  OEPTHS  OF  THE  SIX 

C  SELECTED  ISOTHERMS,  AS  CHOSEN  IN  SUBROUTINE  INTERP 

C 

C  KN — THE  SUBSCRIPT  OF  THE  DEPTH  CHOSEN  AS  THE  INITIAL 

C  APPROXIMATION  OF  THE  CENTER  OF  THE  THERMOCLINE 

C 

C  ZT — THE  DEPTH  OF  THE  TOP  OF  THE  THERMOCLINE 

C 

C  ZC — THE  DEPTH  OF  THE  CENTER  OF  THE  THERMOCLINE 

C 

C  ZBT — THE  DEPTH  OF  THE  BOTTOM  OF  THE  THERMOCLINE 

C 

C  J — USED  AS  THE  SUBSCRIPT  OF  THE  DEPTHS  OBTAINED  FOR 

C  EACH  DATA  CARD 

C 

C  M — THE  GROUPING  NUMBER,  A  MEASURE  OF  THE  NUMBER  OF 

C  INDIVIDUAL  FOURIER  COEFFICIENTS  INCLUDED  IN  EACH 

C  SPECTRAL  BANDWIDTH.   NOTE  THAT  M  MUST  BE  SPECIFIED  IN 

C  THE  MAIN  PROGRAM 

C 

C  AJ,BJ,ZJ — COUNTERS  ASSOCIATED  WITH  THE  CONSECUTIVE 

C  NUMBER  OF  THE  SPECTRAL  DENSITIES,  FOR  USE  IN  GRAPHICAL 

C  OUTPUT 

C 

C  END11,END10,ENDT,ENDC,ETC. — FINAL  SPECTRAL  ESTIMATES 

C  USED  TO  OBTAIN  AN  AVERAGE  POWER  SPECTRUM 

C 

C  SUBROUTINES  CALLED  FROM  MAIN  PROGRAM  ARE 

C  INTERP 

C  NONGA 

C  AVERGE 

C  HANN 

C  SMOOTH 

C  PERIOD 

C  GROUP 

C  RHARM — FOR  DETAILS  CONCERNING  SUBROUTINE  RHARM  SEE  THE 

C  DOCUMENTATION  FOR  THE  IBM  SCIENTIFIC  SUBROUTINE 

C  LIBRARY 


C 


DIMENSION  IT (200) ,TX( 200) ,X( 200) . APER (200 ) ,AJ (200 ) , 
*ZJ(200),BJ(200),Z11(600),Z10(600),Z09(600),Z08(600), 
*Z07(600) ,Z06(600),B11 (600 ) ,B10 (600 ) ,B09( 600) . B08( 600) , 
*BC7( 600), 606(600) ?ZT( 600) ,ZC( 600) ,ZBT( 600) ,AZT(600I, 
*AZC(600) ,AZBT{600),AZIN(600),AZOUT(600),ZIN(600), 
*B0UT(600) ,ARGIN(600) , ARGOUT( 600 ) , INV( 64) , S (64 ) . 
*ENDT(200) , ENDC ( 200 ), ENDB ( 200 ),END1K 200), END1 0(200), 
*END09 ( 200 ) , END08 ( 200 ) , END07 ( 200 ) , END06( 200 ) , T( 200 ) , 
*Z(200) 

REAL  LABL1/1  11  ,/,LABL2/«  10  «/,LABL3/«   9  •/, 
*LABL4/«   8  •/,LABL5/«   7  •/,LABL6/»   6  •/,LABL7/»  TOP* 
*/,LABL8/»MIDL,/,LABL9/«B0TM«/ 

REAL*8  ITI1(12I/"B0X  80D9   AVERAGED  POWER  SPECTRUM  OF 
♦FOURIER     COEFFICIENTS  DATA  INPUT  HAS  MEAN  OF  ZERO  H, 


*S4,ZT»/ 

REAL*8  ITI2! 12)/'80X  80D9   AVERAGED  POWER  SPECTRUM  OF 
♦FOURIER     COEFFICIENTS  OATA  INPUT  HAS  MEAN  OF  ZERO  H, 
*S4,ZCV 

REAL*8  ITI3! 12)/'BOX  80D9   AVERAGED  POWER  SPECTRUM  OF 
♦FOURIER     COEFFICIENTS  DATA  INPUT  HAS  MEAN  OF  ZERO  H, 
*S4,ZBV 

REAL*8  ITITL1  !12)/«BOX  80  FOURIER  ANALYSIS  OF  FLUCTUA 
*TIONS  OF  THE  11  DEGREE  CENTIGRADE  ISOTHERM  RHARM  A 
ANALYSIS'/ 

REAL*8  ITITL2I12)/,B0X  80  FOURIER  ANALYSIS  OF  FLUCTUA 
*TIONS  OF  THE  10  DEGREE  CENTIGRADE  ISOTHERM  RHARM  A 
ANALYSIS'/ 

REAL*8  ITITL3!12)/»B0X  80  FOURIER  ANALYSIS  OF  FLUCTUA 
*TIONS  OF  THE  9  DEGREE  CENTIGRADE  ISOTHERM  RHARM  A 
*NALYSIS»/ 

REAL*8  ITITL4! 12) /'BOX  80  FOURIER  ANALYSIS  OF  FLUCTUA 
♦TIONS  OF  THE  8  DEGREE  CENTIGRADE  ISOTHERM  RHARM  A 
*NALYSISV 

REAL*8  ITITL5( 12) /'BOX  80  FOURIER  ANALYSIS  OF  FLUCTUA 
*TIONS  OF  THE  7  DEGREE  CENTIGRADE  ISOTHERM  RHARM  A 
*NALYSIS*/ 

REAL*8  ITITL6! 12>/»B0X  80  FOURIER  ANALYSIS  OF  FLUCTUA 
♦TIONS  OF  THE  6  DEGREE  CENTIGRADE  ISOTHERM  RHARM  A 
♦  NALYSI  SV 

DO  3  IGROUP=l,2 

DO  1  J=l,512 

READ(5,2)  IITll),  1=1,20) 
!  FORMAT  U3X,20I3) 

CALL  INTERP! I T , 20 , T , Z , ICOUNT, KN, Z 1 1 , Z 10, Z09, Z C8, Z0  7, 
*Z06, J) 

CALL  NONGA(Z,T, ICOUNT , ZT , ZC , ZBT ,KN , J  I 
.  CONTINUE 

CALL  AVERGE  (ZT,AZT) 

CALL    AVERGE    {ZC,AZO 

CALL  AVERGE  (ZBT,AZBT) 

CALL  AVERGE  (Z11,B11) 

CALL  AVERGE  (ZIO.BIO) 

CALL  AVERGE  (Z09.B09) 

CALL  AVERGE  (Z08,B08) 

CALL  AVERGE  (Z07,B07) 

CALL  AVERGE  (Z06,B06) 

CALL  RHARM(AZT,8, INV, S.IFERR) 

CALL  RHARM(AZC,8, INV.S, IFERR) 

CALL  RHARM(AZBT,8,INV,S, IFERR) 

CALL  RHARM(B11,8,  INV, S, IFERR) 

CALL  RHARM<B10,8,  INV, S, IFERR) 

CALL  RHARM(B09,8, INV, S, IFERR) 

CALL  RHARM(B08,8, INV, S, IFERR) 

CALL  RHARM(B07,8,INV,S, IFERR) 

CALL  RHARM(B06,8, INV, S, IFERR) 

CALL  HANN  (AZT.ZT) 

CALL  HANN  <AZC,ZC) 

CALL  HANN( AZBT,ZBT) 

CALL  HANN  <B11,Z11) 

CALL  HANN  <B10,Z10) 

CALL  HANN  (B09.Z09) 

CALL  HANN  (B08,Z08) 

CALL  HANN  (B07,Z07) 

CALL  HANN  (B06,Z06) 

M=A 

CALL  SMOOTH   ( ZT, M , I  TEST , AZT , AJ  ,  I PTS ) 

CALL  SMOOTH   ( ZC , M, IT  EST , AZC, AJ , IPTS ) 

CALL  SMOOTH  ( ZBT , M , I  TEST , AZBT , AJ, IPTS ) 

CALL  SMOOTH(Zll,M,ITEST,Bli,AJ, IPTS) 

CALL  SMOOTH! Z10,M, I  TEST, BIO, A J , IPTS ) 

CALL  SMC0TH(Z09,M,ITEST,B09,AJ,IPTS) 

CALL  SMOOTH! Z0  8, M , I  TEST, B08, A J , IPTS ) 

CALL  SMOOTH! Z07,M, I  TEST , B07, A J , IPTS) 

CALL  SMOOTH! Z06,M, I TEST,B06 , A J . I PTS) 

CALL  PERIODIM,  ITEST, APER, AJ, BJ  ) 

WRITE!6,30) 
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30  FORMAT(lHl) 

WRITE(6,6)  IGROUP 

6  FORMATI//32X,  'SPECTRAL  ESTIMATES  FOR  GROUP  NUMBER  ', 
♦lit//) 

WRITE(6,7) 

7  FORM AT <30X,' PERIOD' ,10X,'TOP' , AX, 'CENTER' ,4X, 
♦•BOTTOM*,/) 

DO  8  1=1, IPTS 

WRITE (6,9)  APER( I ) ,A2T(I) , AZC ( I ) , AZ8T < I ) 
9  FORMAT(30X,F6.2,4X,3F10.2) 

8  CONTINUE 
WRITE(6,30) 
WRITE(6,6)  IGROUP 
WRITE(6,10) 

10  FORMATI30X, 'PERIOD', 7X,'  11  DEG',4X,'10  DEG',5Xt'9  DEC 
♦t/) 

DO  11  I=lfIPTS 

WRITE! 6,9)   APER(  I),B11(I ),B10( I) ,B09(I) 

11  CONTINUE 
WRITE(6,30) 
WRITE(6,6)  IGROUP 
WRITE<6,13) 

13  FORMAT(30X,'PERIOD',8X,'8  DEG',5X,'7  DEG',5X,'6  DEC, 

♦  /) 

DC  14  1=1, IPTS 

14  WRITE<6,9)   APER( I  ),B08(I ),B07( I) ,B06(I) 
DO  26  1=1,  IPTS 

ENDT( < IGROUP-l)*IPTS*I)=AZTU) 

ENDCC  (  IGROUP-  DMPTS+  I  )  =  AZC  (  I  ) 

ENDB( ( IGR0UP-1)^IPTS*I)=AZBT( I) 

END  IK  (IGR0UP-1)*IPTS*I)=B11U) 

END10( UGROUP-l)*IPTS  +  I)  =  B10( I) 

END09( (IGRGUP-1)*IPTS+I)=B09( I) 

END0  8(  (IGROUP-l)*IPTS*I)=B08(  I) 

END07( (IGR0UP-1)MPTS+I)  =  B07( I) 

END06( (IGRCUP-1)*IPTS+I)=B06( I) 
26  CONTINUE 
3  CONTINUE 

CALL  GROUP(AJ,BJ,APER,IPTS,ENDT,ENDC,ENDB,END11,END10, 
♦END09,END08, EN007,END06) 
C      CALL  DRAW  FOR  ISOTHERMS 
C      CALL  DRAW  THERMOCLINES 

C      INCREASE  SIZE  OF  ELEMENTS  OF  ENDT,ENDC,  AND  ENDB 
C      FOR  GRAPHING 

DO  999  1=1, IPTS 

ENDTU  )=1.5+ENDT(  I) 

ENDC(I)=1.5*ENDC( I) 
999  ENDB(I)=ENDB( I )+1.5 

WRITE(6,30) 

CALL  DRAW! IPTS ,B J , END 11 ,0 ,0, LABL1 , IT ITL1 , 1 . ,20.  ,0  ,0  ,  2  , 

♦  2, 5, 8,0, ID 

CALL  DRAW(IPTS,BJ,END10,0,0,LABL2,ITITL2,1.,20.,0,0,2, 

♦  2, 5, 8,0, ID 

CALL  DRAW(IPTS,BJ,END09,0,0,LABL3,ITITL3,1.,20.,0,0,2, 

♦  2, 5, 8,0, ID 

CALL  DRAW( IPTS ,BJ ,END08,0 ,0,LABL4,ITITL*, 1 • ,20. ,0 ,0,2, 

♦  2,5,8,0, ID 

CALL  DRAW( JPTS,BJ,END07,0,0,LABL5,ITITL5,1.,20.,0,0,2, 

♦  2,5,8,0, ID 

CALL  DRAW! IPTS , BJ , END06,0,0,LABL6, IT ITL6, 1  •• 20. ,0,0, 2, 

♦  2,5, 8,0, ID 

CALL  DRAW! IPTS, BJ, ENDT, 0, 0,LABL7, IT  1 1 , 1. , 20. , 0,0, 2, 2, 

♦  5, 8,0, ID 

CALL  DRAWt IPTS ,B J ,ENDC,0, 0,LABL8, I TI2 ,1. ,20. ,0,0,2,2, 

♦  5,8,0, ID 

CALL  DRAWI IPTS ,B J ,ENDB,0, 0, LABL9, IT  13,1., 20. ,0,0,2,2, 

♦  5,8,0, ID 
STOP 

END 


fi? 


SUBROUTINE  INTERPUT,  1PTS  ,  TX  ,  X  ,  I  COUNT  ,  KN,  Z  11  ,  ZIO,  Z09, 
*Z08,ZC7,Z06,KK) 

C  THIS  SUBROUTINE,  CALLED  FROM  THE  MAIN  PROGRAM,  IS  USED 

C  TO  LINEARLY  INTERPOLATE  BETWEEN  OATA  POINTS  GIVEN  AT 

C  FIVF  METER  INCREMENTS  TO  DEPTHS  CORRESPONDING  TO  TEMP. 

C  INCREMENTS  OF  0,1  DEGREE  CENTIGRADE.   THIS  PROGRAM  DOES 

C  NOT  CONSIDER  SMALL  NEAR-SURFACE  TRANSIENTS  OR  INVERSIONS 
C 

C  VARIABLES  USED  IN  THIS  SUBROUTINE 

C  IT — INPUT  TEMPERATURES,  IN  TENTHS  OF  A  DEGREE,  FROM 

C  MAIN  PROGRAM 
C 

C  TX — THE  ARRAY  OF  INTERPOLATED  TEMPERATURES 
C 

C  X--THE  ARRAY  OF  INTERPOLATED  DEPTHS  ASSOCIATED  WITH  TX 
C 

C  ICCUNT — THE  TOTAL  NUMBER  OF  ELEMENTS  IN  X  AND  IN  TX 
C 

C  IPTS — THE  NUMBER  OF  TEMPERATURES  PER  IRACE  AS  READ 

C  FRCM  DATA  CARDS.      IN  THE  CASE  OF  THIS  DATA  IPTS=20 
C 

C  IT1,  IT2,IT3, IT4,IT5,IT6— THE  TEMPERATURES,  IN  TENTHS 

C  OF  A  DEGREE,  OF  THE  ISOTHERMS  WHOSE  DEPTH  FLUCTUATIONS 

C  ARE  OF  INTEREST.   THESE  VALUES  MUST  BE  SPECIFIED  IN 

C  THIS  SUBROUTINE,  AND  ARE  ASSOCIATED  WITH  Z11,Z10,Z09, 

C  Z08,Z07,Z06  RESPECTIVELY  IN  LATER  COMPUTATIONS. 
C 

C  T,TT,ITT — ARRAY  NAMES  ASSIGNED  TO  THE  INTERPOLATED 

C  TEMPERATURES  AT  VARIOUS  POINTS  IN  THE  COMPUTATIONS 
C 

C  Z,IZ — ARRAY  NAMES  ASSIGNED  TO  THE  INTERPOLATED  DEPTHS 
C 

C  ITO — AN  ARRAY  USED  TO  FIND  A  TEMPORARY  TOP,  BOTTOM, 

C  AND  CENTER  OF  THE  THERMOCLINE. 
C 

C  ITTOP — THE  TEMPERATURE  OF  THE  TEMPORARY  TOP  OF  THE 

C  THERMOCLINE.   DEFINED  IN  THIS  SUBPROGRAM  AS  THE  DEPTH 

C  BELOW  WHICH  THE  TEMPERATURE  FIRST  DECREASES  MORE  THAN 

C  0.3  DEGREES  IN  FIVE  METERS  OF  DEPTH. 
C 

C  ITBOT — THE  TEMPERATURE  OF  THE  TEMPORARY  BOTTOM  OF  THE 

C  THERMOCLINE.   DEFINED  AS  THE  FIRST  TEMPERATURE  LESS 

C  TFAN  ITTOP  ABOVE  WHICH  THE  TEMPERATURE  DECREASES  LESS 

C  THAN  0.2  DEGREES  IN  FIVE  METERS  OF  DEPTH 
C 

C  ITC — THE  ARITHMETIC  MEAN  OF  ITTOP  AND  ITBOT,  USED  IN 

C  SLBROUTINE  FINDZC  TO  DETERMINE  THE  SUBSCRIPT  OF  THE 

C  DEPTH  VALUE  CORRESPONDING  TO  THE  ESTIMATED  CENTER  OF 

C  ThE  THERMOCLINE 
C 

C  SUBROUTINES  CALLED  FROM  THIS  SUBROUTINE  ARE 

C  PURGE 

C  FINDZC 

C  FINDZS 


C 


DIMENSION  I T(l)yTX(l) ,X(1) , Zl 1 ( 1 ) , ZIO ( I ) , Z09( 1 ) , Z08 ( 1), 
*Z07( 1) ,Z06(1),T(200) , TTI200) , Z ( 200) , ZZ ( 200 ) , ITT (200), 
*IT0(200) 

ILESS=IPTS-1 

GENERATE  Z(I ) 

DC  999  1=1, IPTS 

Z(  I)  =  5.*( 1-1) 
999  CONTINUE 

IEPS=0 

IT1- 110 

IT2=100 

IT3=90 

ITA=80 
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IT5=70 

IT6=60 

DO    90    II=1,ILESS 

ITCH  II)=IT(IH-1>-IT<II) 

IFUT0(II).LT.-3)    GO    TO    91 

90  CONTINUE 

91  ITTOP=IT(II) 

DO    92    JJ=10tILESS 
IT0(JJ)=IT(JJ*1)-IT( J J) 
IF(IT0(JJ).GT.-2)    GO    TO    93 

92  CONTINUE 

93  ITB0T=IT1JJ) 
ITC=(ITTOP*ITBOT)/2 

C      EXAMINE  ITU)  VERSUS  ITC  1*11  •   IF  SHALLOWER  TEMP, 
C      IS  GREATER  THAN  NEXT  TEMPERATURE  BELOW  IT  IN  THE 
C      WATER  COLUMNt  INTERPOLATE.   IF  NOT,  STORE, 
DO  1000  I=1,IPTS 
J=I+IEPS 
K=0 

IF(I.EQ.IPTS)  GO  TO  71 
IFUTU).GT.IT(I  +  l))  GO  TO  5 
ITTIJI-ITI II 
ZZ(J)=Z(I) 
GO  TO  1000 
5  ITDIF=IT(I  l-ITU  +  ll 
ZDIF=5./ITDIF 
MAXMTDIF  +  J-1 
DO  1001  JJ=J,MAX 
ZZ< J J )=( J J-J)*ZDIF+Z< I) 
ITT(JJ)  =  ITU)-K 
K=K*1 
1001  CONTINUE 

IEPS=MAX-I 
1000  CONTINUE 

INT=ILESS+J 
71  ITT<J)=IT(I) 
ZZ(J)=Z(I) 
100  DO  8  1=1, J 
TTIII-ITTf II 
8  CONTINUE 
ICOUNT=J-l 

CALL  PURGE  (TT ,ZZ. ICOUNT, TX, X ) 
CALL  FINDZC(TX,ITC,KN,ICOUNT) 

CALL  FINDZSCTX.X,  IT1,  IT2,  IT3, IT A, IT5, IT6f Z 11 ,Z10, Z09, 
*Z08,ZC7,Z06,KKj 
20  RETURN 
END 


SUBROUTINE  PURGE  ( TT,ZZ , ICOUNT,TX,X) 

C  SUBROUTINE  PURGE  TAKES  AS  INPUT  THE  ARRAY  OF  TEMPERATURE 

C  AND  DEPTH  INFORMATION  FROM  SUBROUTINE  INTERP  AND 

C  EXAMINES  IT  FOR  DUPLICATION  OF  TEMPERATURE,  CONSIDERING 

C  A  POSITIVE  TEMPERATURE-DEPTH  GRADIENT  TO  BE  ISOTHERMAL. 

C  (THIS  WAS  JUSTIFIED  BY  THE  PARTICULAR  DATA  INVOLVED  IN 

C  THE  WRITING  OF  THIS  SUBPROGRAM.)   THE  SUBROUTINE  STORES 

C  EACH  TEMPERATURE  IN  THE  TRACE  AT  THE  SHALLOWEST  DEPTH 

C  AT  WHICH  IT  APPEARS. 
C 

C  VARIABLES  USED  IN  THIS  SUBPROGRAM 

C  TT,TX,ZZ,X— AS  DEFINED  IN  SUBROUTINE  INTERP 
C 

C  ICOUNT — AS  OUTPUT,  THE  NUMBER  OF  TEMPERATURES  OR  DEPTH 

C  PAIRS  REMAINING  AFTER  DUPLICATIONS  HAVE  BEEN  REMOVED. 


C 


DIMENSION  TT(l),TX(l),ZZ(l),xm 

TX(1)=TT(1) 

X(1)  =  ZZU) 

J=2 
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DO  1  I=l,ICOUNT 
IF(TT(I*i).GE«TT( It)  GO  TO  1 
TX(J)=TT(I+1) 
X(J)=ZZ(I*1) 

J  =  J*1 
CONTINUE 
ICOUNT  =  J-2 
RETURN 
END 


SUBROUTINE  FINDZS(TX,X, IT1,IT2,IT3,IT4,IT5,IT6,Z11, 
*Z10.Z09,Z08,Z07?Z06,KK) 

C  THIS  SUBROUTINE  TAKES  THE  PURGED  ARRAYS  OF  TEMPERATURES 

C  AND  DEPTHS  AND  FINDS  THE  DEPTHS  OF  ISOTHERMS  SELECTED 

C  IN  SUBROUTINE  INTERP 
C 

C  VARIABLES  USED  IN  THIS  SUBROUTINE 

C  TX,X — AS  DEFINED  IN  SUBROUTINE  INTERP 
C 

C  IT1.IT2, IT3,IT4,IT5.IT6 — AS  DEFINED  IN  SUBROUTINE 

C  INTERP 
C 

C  Z11,Z10»Z09,Z08,Z07,Z06 — AS  DEFINED  IN  MAIN  PROGRAM 
C 

C  KK — THE  SUBSCRIPT  OF  THE  DEPTH  VALUE  IN  THE  2**N  ARRAY 


C 


DIMENSION  TX<l),xm  ,Z1U  L)VZ10(1  ),Z09(1)  ,Z08(  I), 
*Z07(1),Z06<1) 
DO  1  1=1,100 
IFITX(I).LE.ITl)  GO  TO  2 

1  CONTINUE 

2  Z11(KK)=X( I) 
DO  3  J»I«100 
IF(TX(J).LE.IT2)  GO  TO  4 

3  CONTINUE 

4  Z10(KK)=X(J) 
DO  5  K=J,100 
IF(TXCK).LE.IT3)  GO  TO  6 

5  CONTINUE 

6  Z09(KK)=X<K) 
DO  7  L=Kf100 
IF(TX(L).LE.IT4)  GO  TO  8 

7  CONTINUE 

8  Z08(KK)=X(L) 
DO  9  M=L,100 
IF(TX<M).LE.IT5)  GO  TO  10 

9  CONTINUE 

10  Z07(KK)=XIM) 
DC  11  MM=M,100 
IF(TX(MM).LE.IT6)  GO  TO  12 

11  CONTINUE 

12  Z06(KK)=X«MM) 
RETURN 

END 


SUBROUTINE  FINDZC (TX, ITCtKN. ICOUNT) 

C  THIS  SUBROUTINE  DETERMINES  THE  SUBSCRIPT  OF  THE  DEPTH 

C  CORRESPONDING  TO  ITC  AS  FOUND  IN  SUBROUTINE  INTERP 
C 

C  VARIABLES  USED  IN  THIS  SUBROUTINE  ARE 

C  TX.ITC— AS  DEFINED  IN  SUBROUTINE  INTERP 
C 

C  ICOUNT — AS  DEFINED  IN  SUBROUTINE  PURGE 


C 


AC 


C      KN— THE  SUBSCRIPT  CF  THE  DEPTH  ASSOCIATED  WITH  ITC 
C 

DIMENSION  TX(1  ) 
DO  1  I=l,ICOUNT 
IF(TX( I) .LE.ITC)  GO  TO  11 
1  CONTINUE 
11  KN=I 
RETURN 
END 


SUBROUTINE  NONGA( Z.T, ICOUNT,ZT, ZC ,ZBT ,KN, JL) 

C  THIS  SUBROUTINE  TAKES  THE  INPUT  ARRAYS  OF  TEMPERATURE 

C  AND  DEPTH,  DIVIDES  THE  DEPTHS  AT  THE  DEPTH  CORRESPONDING 

C  TO  ITC  AND  PRODUCES  TWO  SYMMETRICAL  DEPTH  ARRAYS,  ONE  A 

C  'REFLECTION*  OF  THE  DEPTHS  FROM  Z<KN)  TO  THE  SURFACE  AND 

C  THE  OTHER  A  'REFLECTION*  OF  DEPTHS  FROM  Z(KN)  TO  THE 

C  BOTTOM  CF  THE  TRACE. 
C 

C  VARIABLES  USED  IN  THIS  SUBROUTINE  ARE 

C  Z,T— AS  DEFINED  IN  SUBROUTINE  INTERP 
C 

C  ICOUNT — AS  DEFINED  IN  SUBROUTINE  PURGE 
C 

C  ZT.ZCZBT — AS  DEFINED  IN  MAIN  PROGRAM 
C 

C  JL— THE  SUBSCRIPT  OF  THE  DEPTHS  ZT.ZCZBT  IN  THE  ARRAY 

C  OF  2**N  TERMS  USED  IN  THE  FOURIER  ANALYSIS 
C 

C  KN — THE  SUBSCRIPT  CORRESPONDING  TO  THE  DEPTH  OF  THE 

C  ESTIMATED  CENTER  OF  THE  THERMOCLINE 
C 

C  Z5 — THE  •REFLECTED'  ARRAY  OF  DEPTHS  FROM  Z(KN)  TO  THE 

C  SURFACE 
C 

C  Z3 — THE  'REFLECTED  ARRAY  OF  DEPTHS  FROM  ZtKN)  TO  THE 

C  BOTTOM  OF  THE  TRACE 
C 

C  NUMBER,KDUMMY,NNN,NUMMER — COUNTERS  USED  IN  THE 

C  PROCESS  OF  PRODUCING  GAUSSIAN  THERMOCLINES  FROM  THE 

C  NON-GAUSSIAN  INPUT  THROUGH  'REFLECTION' 
C 

C  SK — A  WEIGHTING  FACTOR,  WHICH  IS  APPLIED  TO  THE 

C  VARIANCE  TO  COMPENSATE  FOR  A  LACK  OF  ACCURACY  IN 

C  COMPUTATION  DUE  TO  USE  OF  FINITE  DIFFERENCE  TECHNIQUES 

C  RATHER  THAN  INTEGRATION.  ITS  VALUE  IS  DEPENDENT  UPON 

C  THE  TEMPERATURE  INTERVAL  (DELTA  Tl  USED  IN  COMPUTATION 
C 

C  SUBROUTINES  CALLED  FROM  THIS  SUBROUTINE  ARE 

C  GAUSS1 

C  GAUSS2 


C 


DIMENSION  Z(200),Z3(200),Z5(200),T(200),ZC<1) ,ZT(1), 
*ZBT(1) 
IFUT<1)-T(2)).GE.10.)  SK=  1.30000 
IF((T{1)-TI2M.LE.5.)  SK=  1.470000 
NUMBER=ICOUNT  -KN+1 
KDUMMY=ICOUNT-KN 
NNN=(2*( ICOUNT-KN)  )  +  l 
ZC(JL»=Z(KM 
DO  73  11=1, KN 
Z5U1)  =  ZU1) 
73  CONTINUE 
IZ1=KN-1 
DO  7  12=1,  IZ1 
X  =  Z<KN)  -  ZCKN-I2) 
Z5CKN+I2)  =  Z(KN)  +  X 
7  CONTINUE 

NUMMER=2*KN-1 


CALL    GAUSSKZ5,NUMMER,IC0UNT,  SK,ZT,ZC,JL) 
4    KKK=KN-1 

DC    74    13=1, NUMBER 

Z3<KDUMMY  ♦  13)  =  Z(KKK*I3) 
74  CONTINUE 

DO  59  I4=1,KDUMMY 

X  =  Z1KN+I4)  -  Z(KN) 

Z3<NUMBER  -  14)  =  Z<KN)  -  X 
59  CONTINUE 

CALL  GAUSS2<Z3,NNN,IC0UNT,SK,ZBT,JL) 

RETURN 

END 


SLBROUTINE  GAUSS  1 (X, NT, ICOUNT.SK, ZT, ZC , JL ) 

C  THIS  SUBROUTINE  FINDS  THE  DEPTHS  OF  THE  TOP  ANO  CENTER 

C  OF  THE  THERMOCLINE  THROUGH  GAUSSIAN  METHODS 
C 

C  VARIABLES  USED  IN  THIS  SUBROUTINE  ARE 

C  X — CORRESPONDS  TO  THE  Z5  ARRAY  OF  SUBROUTINE  NONGA 
C 

C  NT — THE  NUMBER  OF  ELEMENTS  IN  THE  ARRAY  OF  X 
C 

C  Zl — THE  DEPTHS  USED  IN  COMPUTATION  OF  THE  GAUSSIAN 

C  THERMOCLINE.   BY  DEFINITION,  THE  AVERAGE  VALUE  OF 

C  CONSECUTIVE  DEPTHS  IN  THE  X  ARRAY,  THE  MIDPOINT 

C  OF  »DELTA  T» ,  AT  WHICH  THE  DEPTHS  MAKE  THEIR 

C  CONTRIBUTION 
C 

C  12 — THE  ARRAY  OF  VALUES  CORRESPONDING  TO  THE  SQUARES 

C  OF  THE  VALUES  OF  Zl 
C 

C  JL — AS  DEFINED  IN  NONGA 
C 

C  SK — AS  DEFINEO  IN  NONGA 
C 

C  SIGSQR — THE  VARIANCE  CF  THE  ARRAY  OF  Zl 
C 

C  SUM1 — THE  SUM  OF  THE  ARRAY  OF  Zl 
C 

C  SUM2 — THE  SUM  OF  THE  ARRAY  OF  Z2,  OR  THE  SUM  OF  THE 

C  SQUARES  OF  THE  VALUES  OF  THE  ARRAY  OF  Zl 


C 


DIMENSION  Zl(200),Z2(200),X(i),ZC<l),ZTt 1) 

N-NT-1 

SUM1  =0.0 

SUM2=0.0 

DO  I  I  =  1.N 

ZKI)  =  (XU)*X<I*l))/2. 

Z2(I)=Z1(I)**2 

SUM1=SUM1+Z1(I) 

SUM2=SUM2*Z2U) 

ZC(JL)=SUM1/N 

SIGSQR=f  SUM2/N)-CZCUL)**2) 

ZT(JL)=ZC(JL)-SK*SQRT I SIGSQR) 

RETURN 

END 


SUBROUTINE  GAUSS2 C X, NT, ICOUNT.SK, ZBT, JL) 
C    THIS  SUBROUTINE  FINDS  THE  DEPTH  OF  THE  BOTTOM  OF  THE 
C    THERMOCLINE  THROUGH  GAUSSIAN  METHODS 
C 

C    VARIABLES  USED  IN  THIS  SUBROUTINE  ARE 
C      X — CORRESPONDS  TO  THE  Z3  ARRAY  OF  SUBROUTINE  NONGA 


C 
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C  NT — THE  NUMBER  OF  ELEMENTS  IN  THE  ARRAY  OF  X 

C 

C  JL,SK — AS  DEFINED  IN  SUBROUTINE  NONGA 

C 

C  SUM1,SUM2,Z1,Z2— AS  DEFINED  IN  SUBROUTINE  GAUSS1 

C 

C  SIGSQR— AS  DEFINED  IN  SUBROUTINE  GAUSSi 


C 


DIMENSION  Zl ( 200 ), 12 ( 200 ),ZC( 600 ),X< 1),ZBT<1) 
N=NT-1 
SUM1  =0.0 
SUM2=0.0 
DO  1  I=1,N 

ZKI  )=(X(  I)  +  X(  1  +  1)  )/2. 
Z2(I)=Z1(I)**2 
SUM1  =  SUMH-Z1(I  ) 
1  SUM2=SUM2+Z2(I  ) 
ZC(JL)=SUM1/N 

SIGSQR=(SUM2/N)-(ZCt JL>**2) 
ZBT(JL)=ZC(JL)+SK*SQRT( SIGSQR) 
RETURN 
END 


SUBROUTINE  AVERGE ( AZ I N, AZOUT ) 

C  THIS  SUBROUTINE  TAKES  THE  ARRAY  OF  DEPTHS  OF  A  SELECTED 

C  FEATURE  (ISOTHERM  OR  THERMOCLINE ) t  AND  PRODUCES  AN  ARRAY 

C  WITH  A  MEAN  VALUE  OF  ZERO 
C 

C  VARIABLES  USED  IN  THIS  SUBROUTINE  ARE 

C  AZIN — THE  ARRAY  OF  INPUT  VARIABLES,  WITH  AN  ARBITRARY 

C  MEAN  VALUE 
C 

C  AZOUT — THE  OUTPUT  ARRAY  OF  VARIABLES,  WITH  A  MEAN  OF 

C  ZERO 
C 

C  ZINBAR — THE  SUMMATION  OF  VALUES  OF  AZIN 
C 

C  AZINB — THE  MEAN  OF  VALUES  OF  AZIN 


C 


DIMENSION  AZIN(1),AZ0UT(1) 

ZINBAR=0.0 

DO  1  1=1,512 

1  ZINBAR=ZINBAR+AZIN(I) 
AZINB=ZINBAR/512. 

DO  2  1=1,512 

2  AZOUTU)  =  AZIN(I)-AZINB 
RETURN 

ENO 


SUBROUTINE  HANN   <ZIN,BOUT) 

C  THIS  SUBROUTINE  TAKES  THE  FOURIER  COEFFICIENTS  AS  FOUND 

C  IN  SUBROUTINE  RHARM,  AND  HANNS  THEM,  TO  REDUCE  LEAKAGE 

C  BETWEEN  NEIGHBORING  FREQUENCY  BANDS. 
C 

C  VARIABLES  USED  IN  THIS  SUBROUTINE  ARE 

C  ZIN— THE  ARRAY  OF  FOURIER  COEFFICIENTS  FOUND  BY  RHARM 
C 

C  BOUT— THE  HANNED  FOURIER  COEFFICIENTS 


C 


DIMENSION  ZIN( 1 ) , BOUT ( 1 ) 

DO  1  1=1,2 

B0UT<I)=.5*ZIN<I)-.25*ZIN(I+2) 

1  CONTINUE 

DO  2  1=509,510 
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B0UTU+2)=-.25*ZIN<  I  >*.5*ZIN<  1*21 

CONTINUE 

DO  3  1=3,510 

BCUT(I)=-.25*ZIN< I-2)*.5*ZIN( I)-.25*ZIN< 1*2) 

CONTINUE 
RETURN 

END 


SUBROUTINE  SMOOTH ( ARGIN, M, ITEST.ARGOUT ,ZJ, ICOUNT ) 
C    THIS  SUBROUTINE  TAKES  THE  HANNED  FOURIER  COEFFICIENTS 
C    AND  AVERAGES  THEM  OVER  FREQUENCY  BANDS  OF  SPECIFIED 
C    WIDTH.   THE  NUMBER  OF  FREQUENCIES  INCLUDED  IN  EACH  BAND 
C    IS  EQUAL  TO  M,  AS  REAO  IN  BY  THE  MAIN  PROGRAM  BEFORE 
C    THIS  SUBROUTINE  IS  CALLED. 
C 

C    VARIABLES  USED  IN  THIS  SUBROUTINE  ARE 
C      ARGIN — THE  ARRAY  OF  HANNED  FOURIER  COEFFICIENTS 
C 

C      M— A  GROUPING  FACTOR,  AS  DISCUSSED  ABOVE 
C 

C      ITEST — THE  NUMBER  OF  GROUPS  OF  2M  FOURIER  COEFFICIENTS 
C      WHICH  IS  ALSO  ONE  LESS  THAN  THE  TOTAL  NUMBER  OF 
C      SMOOTHED  FREQUENCY  BANDS 
C 

C      ARGOUT — THE  ARRAY  OF  ITEST^l  FREQUENCY  BANDS 
C 

C      ZJ — AS  DEFINED  IN  MAIN  PROGRAM 
C 

C      ICOUNT  — THE  TOTAL  NUMBER  OF  FREQUENCY  BANDS  IN  ARGOUT 
C 

C      FACTOR — A  WEIGHTING  FACTOR  USED  IN  SMOOTHING 
C 

C      N — THE  NUMBER  OF  VALUES  IN  THE  ARRAY  OF  DEPTHS 
C      SUBJECTED  TO  THE  FOURIER  ANALYSIS 
C 

DIMENSION  ARGIN(l), ARGOUT! 1),ZJ( 1) 

CM=M 

N=512 

FACTOR=N/CM 

ITEST=N/(2*M) 

JTEST=(N-1)/M 

IFUTEST*2.GE.JTEST)    GO    TO    3 

ITEST=ITESm 
3    CONTINUE 

IC0UNT=ITEST+1 

D01  KK=1,512 

ARGOUT (KK)=0.0 

ARGIN(KK)=ARGIN(KK)**2 

1  CONTINUE 
PART=0.0 

DO  2  JJ=3,M 

2  PART  =  PART«-ARGIN(JJ) 

ARGOUT ( 1 ) =FACTOR* ( ARG I N II) ♦<  2 . *P  ART ) ) 

ZJ<1>*1. 

DO  100  11=2, ITEST 

MM=( ( I2*II)-3)*M)*2 

MST0P=U2*II1-11*M*1 

DO    101    IJ=MM,MST0P 
101    ARGOUTUI)=ARGOUT(II)+ARGIN(IJ) 

ARGOUTIII )=FACTOR*ARGOUT(II) 
100  ZJUI)=II 

MN=MM*M 

ZJ(ICOUNT)=ICOUNT 

IDIFF=511-((ITEST*2)-1I*M 

IF(2*(IDIFF/2).EQ.IDIFF)    GO    TO    30 

IDIFF=IDIFF*1 
30  PRDIFF=IDIFF/2 

DIFF=512./PRDIFF 


DO  102  JK=MN,512 
102  ARGOUT(ICOUNT)=ARGOUT(ICOUNT)*ARGIN(JK) 
ARGOUT(ICOUNT)=DIFF*ARG0UT(IC0UNT> 
RETURN 
END 


SUBROUTINE  PERIOD(M, I  TEST , APER , A J ,BJ ) 

C  THIS  SUBROUTINE  IS  USED  TO  GENERATE  THE  PERIOOS 

C  ASSOCIATED  WITH  THE  SMOOTHED  SPECTRAL  BANDWIDTHS.  ANO  TO 

C  GENERATE  THE  FINAL  SET  OF  COUNTERS  FOR  USE  IN  GRAPHICAL 

C  DISPLAY 
C 

C  VARIABLES  USED  IN  THIS  SUBROUTINE  ARE 

C  M — AS  DEFINED  IN  SUBROUTINE  SMOOTH 
C 

C  ITEST — AS  DEFINED  IN  SUBROUTINE  SMOOTH 
C 

C  APER — THE  ARRAY  OF  PERIODS  ASSOCIATED  WITH  THE  ARRAY 

C  OF  SMOOTHEC  BANDWIDTHS 
C 

C  AJ — AS  DEFINED  IN  MAIN  PROGRAM 
C 

C  BJ — DEFINED  AS  THE  NATURAL  LOGARITHM  OF  AJ  FOR  USE  IN 

C  PRODUCING  A  SEMILOG  GRAPHICAL  OUTPUT 


C 


DIMENSION  APER(l), AJ(1),BJ(1) 

IPTS=ITEST+1 

APER(1)=(2.*512.)/M 

DO    1    1=2, ITEST 

APER(I)=512./((I-1)*M) 

APER ( I PTS ) =2048. /( 1024. -( 512. -(  ( 2* ITEST* 1)*M) )) 

DO  3  1=1, IPTS 

AJ(I  )=512./(M*APER(I) ) 

DC  2  1=1, IPTS 

BJ(I )=ALOG(AJ(I) ) 

RETURN 

END 


SUBROUTINE  GROUP ( AJ,B J , APER, I PTS, ENDT ,ENDC,ENDB,  END11 , 
♦END10, END09, EN008 , END07, END06 ) 

C  THIS  SUBROUTINE  IS  USED  TO  AVERAGE  THE  RAW  SPECTRA 

C  PRODUCED  IN  •IGROUP'  TIMES  THROUGH  THE  MAIN  PROGRAM  , 

C  AND  TO  PRINT  THE  VALUES  ASSOCIATED  WITH  THESE  AVERAGED 

C  SPECTRA 
C 

C  VARIABLES  USED  IN  THIS  SUBROUTINE  ARE 

C  AJ,BJ,APER — AS  DEFINED  IN  SUBROUTINE  PERIOD 
C 

C  IPTS — THE  TOTAL  NUMBER  OF  SPECTRAL  BANDWIDTHS  COMPUTED 

C  IN  SUBROUTINE  SMOOTH 
C 

C  END11.END10, ENDT, ENDC, ETC. — AS  DEFINED  IN  MAIN  PROGRAM 


C 


DIMENSION  AJ(1 ),BJ(1) ,APER(1»,ENDT( I)  , ENDC ( 1 ) , ENDB {  1  )  , 
♦END11C 1), END 10(1) , END09 ( 1 ) , END0  8< 1),END07( 1),END06( I) 
DO  1  1=1, IPTS 

ENOTd  )=.5*(ENDT(  I  t  +  ENDTt  I+IPTS  I ) 
ENDC  <I)  =  .5*<  ENDC  (  I  )*ENOCC  I  +  IPTSM 
ENDBI I l=.5*( ENDB( I )+ENDB( I+IPTS) ) 
END1K  I)=.5*(ENDlim+ENDll(I*IPTSn 
END10(  I)  =  .5*(EN010m*END10(  I  +  IPTSM 
END09I  I)=.5*<END09U)+END09(  I  +  IPTSM 
E  ND08 ( I ) = . 5* ( ENOO  8(1) +END08 ( I ♦ I PT  S ) ) 
END07I I)=.5*(END07(I)+END07(I+IPTSM 
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END06(  I)=.5*(END0  6( I )*END06( I+IPTS)) 
1    CONTINUE 
WRITE(6, 30) 
30    FORMATdHl  ) 
WRITE(6,3) 

3  FORMAT(/M0X, 'AVERAGEO    SPECTRAL    ESTIMATES1,/) 
WRITE(6,4) 

4  FCRMAT(3CX,  •  NR  •  ,4X,» PERIOD*  .lOX^TOP*  ,4X  ,  »  CENTER  •  ,  AX  , 
♦•BOTTOM*, /) 

DO    5    I=1,IPTS 

5  WRITE(6,6)     AJ(  I),  APERU  )  ,ENDT(I  )  ,ENDC(  I)  ,ENOB(  I) 

6  F0RMAT(28X.F4.1,4X,F6.2,4X,3F10.2) 
WRITE(6,30) 

WRITE(6,3) 
WRITE<6,7) 

7  FORMATOOXt'NR*  ,4X,  'PERIOD',  7X,«11    DEG,,4X,,10    DEG», 
*5X,»9    DEG»./) 

DO    8    I=1,IPTS 

8  WRITE(6,6)    AJ(  I)  , APERU)  ,END11(I)  ,EN010(  I),END09(  I) 
WRITE(6,30) 

WRITE<6,3) 
WRITF(6,10) 
10    F0RMAT(30X, •  NR  • , AX, 'PERIOD", 8X, «8    DEG»,5X,«7    DEG*,5X, 
*»6    DEG« ,/) 
DO    9    I=1,IPTS 

9  WRITE<6,6)    AJ( I),APER( I ),END08( I),EN007( I ),END06( I) 
RETURN 

END 
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